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abstract  (Maximum  200  words) 

In  this  research,  we  have  attempted  to  characterize  the  anisotropy  in 
particle  arrangement  and  particle  interactions  that  develop  when  a  granular 
material  is  deformed  and  to  determine  the  influence  of  this  anisotropy  on  the 
subsequent  behavior  of  the  material.  The  research  has  dealt  with  the  behavior 
of  an  idealized  granular  material  consisting  of  glass  spheres  and  has  involved 
the  close  interaction  between  physical  experiments,  numerical  simulations,  and  the 
development  of  theory. 

We  have  focused  on:  (1)  Understanding  the  relationship  between  the 
developing  anisotropy  of  the  material  and  the  changes  in  volume  observed  in 
aeformations  in  which  the  mean  stress  is  held  constant.  In  particular,  we  have 
■onsidered  the  volume  change  of  a  hollow  cylindrical  sample  in  triaxial 
compression  and  extension;  (2)  Using  elastic  waves  to  probe  the  evolution  of  the 
internal  structure  of  the  granular  material  as  it  is  being  deformed;  and  (3) 
Developing  constitutive  relations  for  general  stress  paths  in  the  triaxial/ 
torsional  device  through  the  combined  efforts  of  theoretical  modeling  and 
numerical  simulation. 


SUBJECT  TERMS 

Granular  material,  fabric,  volume  change,  boundary  effects, 
wave  propagation,  anisotropic  elasticity,  numerical  simulatio 
"rlaxial  compression,  polycrystalline  plasticity 


SECURITY  CLASSIFICATION 
OF  REPORT 

ICIASSIFIED 


18.  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 

UNCLASSIFIED 


19. 


security  classification 
OF  abstract 


UNCLASSIFIED 


15.  NUMBER  OF  PAGES 

107 


tS.  PRICE  CODE 


20.  LIMITATION  OF  ABSTRAa 

SAR 


Standard  ‘omi  298  (Rev  2-89' 

o-c  D.  Sic  ‘5 

:-it  _ 


754G-0-.-280-5500 


Table  of  Contents 


Introduction 

J.  T.  Jenkins . ii 

Publications  and  Theses . vii 

Participants . x 

1 .  Influence  of  the  Specimen  Boundary  and  Shape  on  Volumetric  Be¬ 
havior  of  Granular  Materials 

I.  Ishibashi  and  J.  W.  Choi . 1 

2.  Volume  Change  in  Inhomogeneous  Samples 

P.  A.  Cundall . 16 

3.  Anisotropic  Elastic  Constants  of  a  Granular  Assembly  from  Wave 
Velocity  Measurements 

I.  Ishibashi  and  T.  K.  Agarwal . 25 


4.  Determination  of  the  Principal  Acoustic  Axes  in  Anisotropic  Solids 
from  Wavespeed  Measurements 

B.  Castagnede  and  J.  T.  Jenkins . 38 

5.  Anisotropic  Elasticity  for  Random  Arrays  of  Identical  Spheres 

J.  T.  Jenkins . 48 

6.  Evolution  of  Elastic  Moduli  in  a  Deforming  Granular  Assembly 

P.  A.  Cundall . 59 

7.  Measurement  of  Local  Strainrate  in  a  TRUBAL  Sample 

P.  A.  Cundall . 67 

8.  Mean-field  Inelastic  Behavior  of  Random  Arrays  of  Identical 
Spheres 

J.  T.  Jenkins  and  O.  D.  L.  Strack . 78 

9.  The  Plastic  Anisotropy  and  Average  Plastic  Spin  of  Planar  Poly¬ 
crystalline  Aggregates 

V.  Prantil,  J.  T.  Jenkins,  and  P.  R.  Dawson . 94 


'es 


□  oisiL 


Introduction 

James  T.  Jenkins 


This  final  report  provides  an  overview  of  the  research  accomplished  dmring  the 
course  of  the  grant.  It  also  serves  as  access  to  the  publications  resulting  firom  the 
grant.  In  these  publications,  all  aspects  of  the  activity  are  elaborated  upon.  In  this 
Introduction,  we  provide  a  summary  of  the  accomplishments  and  indicate  where 
more  information  can  be  found  in  the  Chapters  that  follow. 

In  the  research,  we  have  attempted  to  characterize  the  anisotropy  in  particle 
arrangement  and  particle  interactions  that  develop  when  a  granular  material  is  de¬ 
formed  and  to  determine  the  influence  of  this  anisotropy  on  the  subsequent  behavior 
of  the  material.  The  research  has  dealt  with  the  behavior  of  an  idealized  granular 
material  consisting  of  glass  spheres  and  has  involved  the  close  interaction  between 
physical  experiments,  niunerical  simulations,  and  the  development  of  theory. 

We  have  focused  on:  (1)  Understanding  the  relationship  between  the  develop¬ 
ing  anisotropy  of  the  material  and  the  changes  in  volume  observed  in  deformations 
in  which  the  mean  stress  is  held  constant.  (2)  Using  elastic  waves  to  probe  the  evo¬ 
lution  of  the  internal  structure  of  the  granular  material  as  it  is  being  deformed;  and 
(3)  Developing  constitutive  relations  for  general  stress  paths  in  the  triaxial/torsional 
device  through  the  combined  efforts  of  theoretical  modeling  and  numerical  simula¬ 
tion. 


Volume  Change 

The  goal  of  the  work  on  volume  change  was  to  develop  an  understanding  of 
how  various  types  of  anisotropy  induced  in  a  granular  material  influence  its  vol¬ 
ume  change  in  subsequent  deformations.  Such  an  understanding  should  lead  to  a 
quantitative  measure  of  the  liquefaction  potential  of  a  saturated  material. 

Much  of  our  effort  has  been  to  understand  the  volume  change  of  a  hollow  cylin¬ 
drical  sample  held  at  constant  pressure  in  a  triaxial/torsion  device.  The  sample  was 
found  to  expand  when  subjected  to  triaxial  compression  and  contract  when  sub¬ 
jected  to  triaxial  extension.  The  theory  predicted  such  behavior  only  if  there  was 
anisotropy  in  the  orientational  distribution  of  contacts.  Consequently,  an  effort  was 
made  in  the  simulations  to  establish  an  initial  anisotropy  in  the  orientational  dis¬ 
tribution  of  contacts  alone.  However,  any  such  initial  anisotropy  disappeared  when 
the  material  was  loaded  with  an  isotropic  stress  and  the  particles  rearranged  to  ac¬ 
commodate  it.  The  indication  was  that  depositional  anisotropy  was  not  responsible 
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for  the  difference  in  the  volume  change  in  triaxial  compression  and  extension. 

In  Chapter  One,  Ishibashi  and  Choi  describe  how  it  was  determined  that  the 
methods  used  to  prepare  the  sample  might  have  induced  an  anisotropy  that  was 
impossible  to  reproduce  in  the  numerical  simulations.  During  the  construction  of 
the  sample,  regions  of  ordered  grains  may  form  near  the  boundaries  of  the  mold 
resulting  in  a  strong  but  localized  anisotropy.  Also,  and  perhaps  more  importantly, 
because  of  the  difference  between  the  boundaries  at  the  ends  and  sides  of  the  hollow 
cylindrical  sample,  it  is  likely  that  triaxial  compression  is  strain  controlled  while  tri¬ 
axial  extension  is  stress  controlled.  Consequently,  there  is  an  anisotropy  associated 
with  the  experimental  device. 

Strain  control  in  triaxial  compression  occurs  because  the  rigid  platens  at  top 
and  bottom  can  exert  contact  forces  of  essentially  arbitrary  magnitude  to  impose 
the  specified  displacement.  As  a  consequence,  they  support  the  development  of 
load  bearing  particle  chains  along  the  axis  of  the  sample.  In  contrast,  the  lat¬ 
eral  membranes  are  limited  in  the  magnitude  of  the  contact  force  that  they  can 
support.  The  particle  displacements  must  accommodate  themselves  to  this,  and, 
consequently,  they  are  not  controlled.  In  triaxial  extension,  load  bearing  particle 
chains  are  suppressed  near  the  boundary.  Expansion  of  the  sample  is  assumed  to 
be  associated  with  anisotropy  aligned  in  the  major  principal  direction.  In  triaxial 
compression,  this  anisotropy  is  allowed  to  develop;  in  triaxial  extension,  it  is  not. 

In  Chapter  Two,  Cundall  employs  the  numerical  program  FLAG  to  investigate 
the  response  of  an  isotropic,  linear  elastic  material  with  elastic  moduli  that  fluctuate 
in  space  about  their  average  values.  Holding  the  mean  stress  constant  and  using 
combinations  of  stress  controlled  and  strain  controlled  boundaries,  he  determines 
the  possible  deformations  of  this  isotropic  but  statistically  inhomogeneous  material. 
He  finds  that  volume  changes  inevitably  occur. 

The  implications  for  soils  testing  of  hollow  samples  confined  by  lateral  mem¬ 
branes  and  rigid  end  platens  are  rather  profound.  For  these,  stress  paths  between 
triaxial  compression  and  extension  involve  a  mixture  of  stress  control  and  strain  con¬ 
trol.  Consequently,  along  such  paths  information  about  the  confinement  is  being 
confused  with  information  about  the  material  behavior.  The  solution  is  to  employ 
thick-walled  samples  or  to  provide  additional  stiffening  to  the  lateral  boimdaxies. 
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Wave  Propagation 


The  research  on  wave  propagation  has  as  it  goal  the  determination  of  the 
anisotropy  of  the  granular  material  from  the  measurements  of  the  speeds  of  shear 
waves  and  pressure  waves  propagating  through  the  material.  This  should  lead  to 
the  capabihty  of  predicting  the  anisotropy  of  a  granular  material  in  the  field. 

Experiments  involving  wave  propagation  through  a  granular  material  during 
its  deformation  were  carried  out  in  a  true  triaxial/torsion  device  equipped  as  a 
resonzmt  column  and  in  a  cubical  cell  with  bending  bimorphs  as  transducers.  In 
both  configurations  the  development  of  anisotropy  with  deformation  could  be  de¬ 
termined  from  the  wave  speed  measurements.  In  the  cubical  cell,  the  evolution  of 
the  complete  set  of  elastic  moduli  could  be  followed.  This  research  is  described  in 
detail  in  Chapter  Three.  The  interpretation  of  the  experiments  was  facilitated  by 
the  simultaneous  development  of  theory  and  numerical  simulations. 

As  outlined  in  Chapter  Four,  algorithms  were  developed  for  the  optimal  de¬ 
termination  of  the  elastic  constants  of  an  anisotropic  material  from  wave  speed 
measurements.  These  are  most  easily  implemented  when  the  axes  of  anisotropy  of 
the  material  is  known  beforehand,  but  can  be  applied  to  determine  these  axes  in 
materials  of  the  common  symmetry  classes.  In  the  experiments  on  granular  mate¬ 
rials  carried  out  so  far,  the  algorithms  have  been  employed  to  determine  the  elastic 
constants  of  a  material  in  which  a  transverse  isotropy  is  developing. 

The  theoretical  calculations  of  Chapter  Five,  based  on  the  mean  field  assump¬ 
tion,  were  used  to  determine  the  elzistic  constants  of  a  random  array  of  spheres  in 
which  the  distribution  of  contacts  exhibited  transverse  isotropy.  The  mean  field 
assumption  usually  equates  the  rotation  of  particles  and  the  displacement  of  the 
centers  of  neighboring  particles  to  the  average  strain  and  rotation  of  the  sample. 
However,  the  calculations  show  that  when  the  contact  distribution  is  anisotropic,  in 
order  to  maintain  a  symmetric  stress,  it  is  necessary  to  permit  the  average  rotations 
of  the  particles  to  differ  from  the  average  rotation  of  the  sample.  Using  the  resulting 
theory,  wave  speeds  could  be  calculated  and  compared  with  numerical  simulations 
and  experiment. 

The  numerical  simulation  program  TRUBAL  was  used  to  provide  the  effec¬ 
tive  moduli  governing  the  wave  propagation  in  an  isotropic  array.  As  discussed  in 
Chapter  Six,  these  numericjil  simulations  were  in  good  agreement  with  the  mod¬ 
uli  measured  in  the  corresponding  experiments,  but  the  shear  modulus  was  about 
one  third  that  predicted  by  theory.  In  an  effort  to  understand  this  discrepancy, 
the  distribution  of  the  normal  component  of  the  contact  force  was  obtained  from 
the  simulation  and  found  to  be  exponential.  However  the  difference  between  the 


measured  and  predicted  wave  speed  did  not  have  its  origin  in  this  distribution. 
What  did  become  clear  during  the  course  of  the  simulation  was  the  extent  to  which 
the  detailed  motions  of  individual  particles  differed  from  those  consistent  with  the 
mean-field  assumption. 

Chapter  Seven  describes  an  attempt  to  characterize  the  local  strains  and  rota¬ 
tions  and  to  compare  them  to  their  average  values.  TRUBAL  is  used  together  with 
a  least-squares  procedure  to  determine  the  local  strain,  rotation,  and  displacement 
of  the  particles  in  a  simulation.  Increments  in  the  average  strain  and  rotation  ten¬ 
sors  calculated  from  the  local  values  are  compared  with  increments  in  the  overall 
strain  and  rotation  tensors.  There  is  a  surprisingly  large  difference  between  two 
kinds  of  averages. 

We  are  presently  carrying  out  similar  simulations  on  larger  samples  in  an  effort 
to  understand  this  discrepancy.  In  order  to  make  correct  quantitatively  predictions 
of  wave  speeds  as  the  anisotropy  of  a  deforming  granular  material  evolves,  it  seems 
necessary  to  characterize  the  deviations  of  the  local  particle  displacements  and  rota¬ 
tions  from  those  associated  with  the  average  strain  and  rotation  and  to  incorporate 
measures  of  these  deviations  into  the  theory. 

Constitutive  Relations 

The  mean-field  assumption  was  also  used  to  determine  the  response  of  the  ma¬ 
terial  when  is  was  subjected  to  triaxial  compression  and  extension.  A  simplification 
of  the  Hertz-Mindlin  contact  law  was  assumed  in  which  the  tangential  component 
of  the  contact  force  was  assumed  to  depend  linearly  and  reversibly  on  the  tangential 
displacement  of  contacting  spheres  until  frictional  slip  occlude.  It  was  assumed  that 
the  shear  strains  were  rather  modest,  never  exceeding  by  much  the  volume  strain 
associated  with  the  isotropic  part  of  the  stress.  In  this  case  it  is  reasonable  to  sup¬ 
pose  that  the  particle  arrangement  remains  unchanged,  although  some  contacts  are 
effectively  lost  as  their  contact  forces  relax  to  zero. 

The  calculation  of  the  stress-strain  relation  for  triaxial  compression  is  given  in 
Chapter  Eight.  The  predicted  response  of  the  material  is  in  qualitative  agreement 
with  the  experiments  and  numerical  simulations,  but  it  is  again  too  stiff.  The  indi¬ 
cation  is  that  here,  too,  measures  of  the  deviations  from  the  mean  field  assumption 
must  be  incorporated  into  the  theory. 

The  theoretical  model  also  permits  the  calculation  of  the  inelastic  parts  of  the 
response  such  as  the  average  plastic  strain,  the  yield  function,  the  plastic  potential, 
and  the  hardening  function,  although  the  utility  of  these  in  describing  the  material 
must  await  our  calculation  of  its  behavior  upon  unloading. 
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An  effort  was  begun  to  formulate  theory  for  shear  strains  large  compared  to 
the  volume  strain  associated  with  the  pressure.  Here  substantial  rearrangement  of 
particle  positions  were  expected  to  accompany  the  particles’  sliding  and  rotation. 
In  an  effort  to  educate  ourselves  to  the  appropriate  type  of  theory  to  employ,  we 
turned  to  the  literature  on  the  polycrystalline  plasticity.  Here  the  elements  of  the 
aggregate  are  single  crystals  that  deform  and  rotate  by  rate-independent  slip.  In 
trying  to  understand  this  subject,  we  found  that  we  could  contribute  to  it. 

Chapter  Nine  describes  a  micromechanical  theory  for  the  development  of  orien¬ 
tational  anisotropy  in  a  simple  planar  polycrystalUne  aggregate  comprised  of  single 
crystals  with  two  shp  planes.  The  introduction  of  the  appropriate  form  of  the 
mean-field  assumption  leads  to  an  equation  of  evolution  for  the  orientation  of  a 
single  crystal  which,  in  turn,  may  be  used  to  determine  the  evolution  of  the  distri¬ 
bution  of  orientations.  Average  properties  can  be  calculated  using  the  orientational 
distribution.  A  quantity  of  particular  interest  is  the  average  plastic  spin.  When  this 
quantity  is  calculated,  its  evolution  is  found  to  depend  on  only  the  second  moment 
of  the  distribution  function,  the  analog  of  the  geometric  fabric  tensor  in  granular 
materials. 

Similar  calculations  have  been  made  in  three  dimensions  for  single  crystals 
with  an  arbitrary  number  of  slip  planes  and  applied  to  face-centered  cubic  crystals. 
The  analogous  theory  for  frictional  slip  of  random  arrays  of  contacting  spheres  is 
presently  being  developed. 
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Chapter  One 


Influence  of  Specimen  Boundary 
and  Shape  on  Volumetric  Behavior 
of  Granular  Materials 

Isao  Ishibashi  James  T.  Jenkins  Jong  W.  Choi 


Abstract 


Laboratory  shear  tests  were  conducted  on  large  and  small  glass 
bead  assemblies  with  both  stress  and  strain  controlled  lateral  bound¬ 
aries  in  both  an  ordinary  solid  cylinder  triaxial  device  and  in  a  hollow 
cylinder  device.  The  specimens  were  sheared  in  triaxial  compression 
and  extension  keeping  the  mean  stress  constant  and  volumetric  strains 
were  monitored.  The  results  showed  that  (1)  the  mold  used  in  prepar¬ 
ing  the  specimen  created  a  strong  but  localized  ordering  of  particles  at 
the  specimen  boundary  and  this  type  of  initial  fabric  had  a  dominant 
effect  on  the  volumetric  behavior  of  the  hollow  cylinder  specimens, 
and  (2)  the  strain  and  the  stress  controlled  boundaries  produced  sig¬ 
nificantly  different  volumetric  strains,  particularly  at  the  initial  stage 
of  the  shearing. 


1  Introduction 

During  the  past  several  years,  the  authors  and  co-researchers  (Chen  et  al.,  1988, 
Ishibashi  et  al.,  1988,  Jenkins,  1988,  Cundall,  1988,  Jenkins  et  al.,  1989,  Ishibashi  et 
al.,  1989,  Ishibashi  and  Agarwal,  1991)  have  been  investigating  anisotropic  behavior 
of  granular  materials,  experimentally,  numerically,  and  theoretically.  A  special  focus 
was  placed  on  the  volumetric  behavior  along  veirious  stress  paths.  The  stress  paths 
that  we  have  had  a  particular  interest  in  are  triaxial  compression  and  extension  tests 


with  constant  mean  stress  (i.e.,  ((7i  +  0-2  +<T3)/3  =  constant).  Typical  experimental 
volumetric  strain  data  on  assemblies  of  glass  spheres  from  hollow  cylinder  tests 
(Chen  et  al.,  1988)  and  ordinary  solid  cylinder  triaxial  tests  (Ishibashi  et  al.,  1989) 
are  shown  in  Figure  1. 


Figure  1:  Volumetric  strains  of  a  fine  glass  bead  assembly  in  triaxial  compression 
and  extension  in  hollow  and  solid  cylinder  specimens  (Chen  et  al.  (1988),  Ishibashi 
et  al.  (1989)). 

Figures  2  and  3  show  computer  simulated  results  using  a  discrete  element  code 
“TRUBAL”  for  nearly  isotropic  and  highly  anisotropic  specimens  (Ishibashi  and 
Agarwal,  1991),  respectively.  Two  questions  were  asked  based  on  these  results;  (1) 
what  is  the  reason  for  the  big  difference  in  volumetric  strain  between  the  triaxial 
compression  and  extension  tests  in  the  hollow  cylinder  specimen?  and  (2)  why  do 
total  dilations  occur  at  the  beginning  of  shear  in  the  triaxial  compression  tests  in  the 
hollow  specimen  and  in  triaxial  extension  in  the  numerical  simulations  but  not  in 
the  solid  cylinder  specimens?  The  researchers  made  two  hypotheses  to  explain  the 
above  phenomena;  (1)  hollow  cylindrical  specimens  have  a  peculiar  inherent  fabric 
due  to  their  relatively  small  thickness,  and  (2)  the  differences  in  the  volumetric 
strain  of  the  specimen  may  be  attributed  to  the  differences  between  the  side  and  end 
boundary  conditions.  The  first  hypothesis  was  based  on  the  observation  that  during 
the  specimen  preparation  an  ordered  arrangement  of  the  particles  around  the  lateral 
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Figure  2:  Computer  simulated  volumetric  strain  for  an  isotropic  specimen  (Ishibashi 
and  Agarwal,  1991). 
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Figure  3:  Computer  simulated  volumetric  strain  for  a  presheared  specimen 
(Ishibashi  and  Agarwal,  1991). 
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boundaries  of  the  molds  might  be  created  and  thus  the  initial  specimen  might  have 
a  peculiar  inherent  fabric.  Because  in  the  hollow  device  the  specimen  is  relatively 
thin  in  terms  of  particle  diameters,  the  eflfect  of  these  ordered  boundary  particles 
might  have  a  significant  influence  on  the  entire  volumetric  behavior.  The  second 
hypothesis  is  based  on  the  fact  that  in  the  triaxial  specimen  the  top  and  bottom 
boundaries  are  made  of  rigid  plates  which  induce  a  uniform  displacement  (strain 
controlled  boundary),  while  the  lateral  confinement  with  thin  rubber  membranes 
imposes  a  stress  that  is  limited  in  magnitude. 

The  authors  reviewed  a  number  of  papers  that  employed  similar  testing  condi¬ 
tions  such  as: 

(a)  Yamada  and  Ishihara,  1979, 

(b)  Haruyama,  1981, 

(c)  Right  et  al.,  1983, 

(d)  Ochiai  and  Lade,  1983, 

(e)  Miura  et  al.,  1986, 

(f)  Sture  et  al.,  1987,  and 

(g)  Lam  and  Tatsuoka,  1988. 

These  testing  devices  included  hollow  cylindrical  specimens  (references  (c)  and  (e) 
above),  solid  cylindrical  specimens  (reference  (g)  above),  and  cubical  specimens 
(references  (a),  (b),  (d)  and  (f)  above).  Those  boundary  conditions  were  stress 
boundaries  (references  (a)  and  (b)  above)  and  mixed  boundaries  (references  (c), 
(d),  (e),  (f)  and  (g)  above).  No  experiments  involving  strain  controlled  boundaries 
were  reported  relevant  to  this  study.  It  should  be  noted,  however,  that  the  numer¬ 
ical  simulations  by  TRUBAL  (Cundall,  1988,  Ishibashi  et  al.,  1989,  Ishibashi  and 
Agarwal,  1991)  are  purely  strain  controlled. 

From  this  review  no  clear  conclusions  could  be  drawn  regarding  the  effect  of 
sample  shape  and  the  boundary  conditions  on  the  volumetric  behavior,  probably 
because  the  importance  of  these  conditions  was  not  recognized  at  the  time  of  the  ex¬ 
periments.  However,  some  interesting  and  relevant  observations  were  made.  A  dense 
hollow  cylinder  specimen  preloaded  in  the  vertical  direction  (Right  et  al.,  1983)  and 
a  loose  specimen  of  glass  beads  in  the  cubical  true  triaxial  device  (Haruyama,  1981) 
showed  a  very  similar  volumetric  behavior  as  that  in  Figure  1.,  while  Yamada  and 
Ishihara’s  cubical  specimens  (1979)  for  a  loose  river  sand  showed  both  contractive 
volume  changes.  Lam  and  Tatsuoka’s  experiment  (1988)  showed  a  cle£ir  effect  of 
stress  and  strain  boundaries  in  volumetric  strain,  although  their  tests  were  not 
under  a  constant  mean  stress. 
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2  Experimental  Procedures 

Two  types  of  glass  bead  mixtures  were  used.  The  first  mixture,  called  the 
“small  bead  assembly”,  was  composed  of  two  sizes  of  glass  beads  with  0.215  mm 
and  0.256  mm  diameters  in  a  number  ratio  of  1  to  1  or  weight  ratio  of  1  to  1.688. 
The  second  mixture,  called  the  “large  bead  assembly” ,  was  composed  of  two  larger 
sizes  of  bead,  3  mm  and  4  mm,  in  a  number  ratio  of  4.728  to  1  or  weight  ratio  of  2 
to  1. 

Two  different  specimen  boundaries  were  created.  For  a  lateral  strain  controlled 
boundary  relatively  rigid  plastic  or  aluminum  liners  were  employed.  The  plastic 
liner  (126.0  mm  high  x  212.9  mm  wide  x  1.143  mm  thick)  is  easily  bent  and  was 
preshaped  to  the  same  curvature  as  the  surface  of  the  specimen.  The  aluminum 
liners  consisted  of  eight  pieces  of  longitudinal  aluminum  plates(126.0  mm  high  x  26.7 
mm  wide  x  1.143  mm  thick)  that  also  had  the  same  curvature  as  the  surface  of  the 
specimen.  The  liners  were  placed  around  the  specimen  just  inside  of  the  membrane 
so  that  the  confining  pressure  would  be  indirectly  imposed  on  the  specimen  via 
the  liners.  The  liners  were  used  only  in  the  solid  specimens  and  not  in  the  hollow 
specimens. 

The  majority  of  the  experiments  were  conducted  in  a  solid  cylindrical  triaxial 
compression  device,  which  accommodated  a  specimen  of  about  71.8  mm  in  diameter 
and  152.4  mm  in  height.  One  series  of  tests  was  performed  in  a  torsional  simple 
shear  apparatus  (Chen  et  al.,  1988),  which  employed  a  specimen  with  a  71.12  mm 
inside  diameter,  101.6  mm  outride  diameter  and  about  142  mm  in  height. 

All  of  the  specimens  were  compacted  (tamped)  in  equal  layers,  four  for  the 
solid  specimens  and  five  for  the  hollow  specimens,  to  an  average  relative  density  of 
60%.  A  burette  was  connected  to  the  bottom  for  measuring  the  volume  change  of 
the  fully  saturated  specimen. 

Each  specimen  was  isotropically  consolidated  from  34.5  kPa  to  138  kPa  prior  to 
being  shewed  along  the  desired  stress  paths.  For  all  tests,  the  mean  effective  stress, 
(<Ti  +  <72  +  <73)73,  was  kept  constant  at  138  kPa  throughout  the  shearing  stage  by 
controlling  the  vertical  stress  and  the  radial  stress  (equal  to  the  tangential  stress). 
The  following  results  on  volumetric  strain  were  obteiined  after  correcting  for  the 
penetration  of  the  membrane  into  the  voids  of  the  glass  beads  at  the  boundary.  A 
detailed  discussion  on  the  membrane  penetration  correction  can  be  found  in  Choi 
and  Ishibashi  (1992). 
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3  Test  Results 


3.1  Vertical  and  Lateral  Strains  Under  Isotropic  Compression 

Figures  4(a)  and  (b)  show  the  relationship  between  the  measured  vertical  strain 
£z  and  the  calculated  radial  strain  Cr  during  isotropic  compressions  fo  the  small  bead 
assemblies  with  lateral  stress  and  strain  controlled  boundaries,  respectively,  in  the 
solid  triaxial  device.  The  radial  strain  €t  was  calculated  from  measured  volumetric 
strain  £„  and  the  vertical  strain  £z  with  an  assumption  of  tangential  strain  £e 
was  equal  to  the  radial  strain  £r.  It  is  generally  accepted  that  this  assumption  is  a 
valid  approximation  for  transversely  anisotropic  media.  These  figures  show  a  nearly 
isotropic  behavior  except  at  the  very  beginning  of  the  compression. 

Figures  5(a)  and  (b)  show  the  same  relationships  as  Figures  4(a)  and  (b)  but  for 
the  large  bead  assemblies.  These  figures  show  a  highly  anisotropic  natme  through¬ 
out  the  entire  compression. 

3.2  Volumetric  Behavior  During  Triaxial  Compression  and  Ex¬ 
tension  Tests 

Figures  6(a)  and  (b)  show  the  relationship  between  volumetric  strain  and  the 
maximum  shear  strain  for  soUd  cylindrical  assemblies  of  small  beads  with  stress  and 
strain  controlled  boundaries  during  triaxial  compression  and  extension  lateral  tests. 
Both  the  triaxial  compression  and  extension  results  are  very  similar. 

Figures  7(a)  and  (b)  show  similar  plots  for  the  large  beads.  In  contrast  to 
Figures  6  (a)  and  (b),  the  extension  test  data  shows  a  large  initial  contraction  while 
the  triaxial  compression  tests  show  very  little  (stress  controlled  lateral  boundary) 
or  no  (strain  controlled  lateral  boundary)  initial  contraction  of  the  speciman.  These 
specimen  seem  to  be  extremely  anisotropic  initially. 

4  Discussion 

4.1  Anisotropy  Due  to  the  Specimen  Preparation  Mold 

The  mold  confinement  diuring  specimen  preparation  could  create  a  strong 
anisotropic  fabric  in  the  vertical  direction  near  the  boundaries  and  the  size  of  this 
region  (the  cross-hatched  zones  in  Figmre  8)  will  increase  as  the  average  grain  size 
increases. 

From  this  point  of  view,  the  large  glass  beads  in  the  solid  triaxial  device  are 
expected  to  provide  a  simileir  volumetric  behavior  to  that  from  the  hollow  cylinder 
specimen  with  the  small  glass  beads.  This  is  because  the  maximum  number  of  the 
large  grains  across  the  diameter  of  the  solid  specimen  is  22  which  is  comparable  to 
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Radial  Strain,  £r  (%) 


Figure  4:  Relationship  between  radial  and  vertical  strains  under  isotropic  compres¬ 
sion  test  for  small  bead  assemblies:  (a)  lateral  stress  boundary,  (b)  lateral  strain 
boundary 


Figure  5:  Relationship  between  radial  and  vertical  strains  under  isotropic  compres¬ 
sion  test  for  large  bead  assemblies:  (a)  lateral  stress  boundary,  (b)  lateral  strain 
boundary.  _ 


70  small  beads  across  the  15.2  mm  thick  wall  of  the  hollow  specimen.  Meanwhile,  in 
the  solid  specimen  the  maximum  number  of  small  galss  beads  across  the  diameter 
is  about  308  so  that  the  mold  confinement  effect  should  be  small.  Figure  9,  which 
shows  volumetric  strains  obtained  from  the  hollow  specimens,  is  consistent  with  the 
above  observation.  That  is.  Figure  9  is  very  similar  to  Figure  7(a). 

As  a  demonstration,  one  solid  specimen  of  small  beads  was  presheared  up  to 
Tmax  =  97  kPa  (or  axjaz  =  3.58)  in  triaxial  compression  and  unloaded  to  the  initial 
isotropic  stress.  A  specimen  with  an  induced  fabric  is  expected  to  exhibit  a  strongly 
anisotropic  volume  change  during  reshearing.  The  results  of  reshearing  are  shown  in 
Figure  10.  The  difference  between  two  orthogonal  stress  paths  (triaxial  compression 
and  extension)  is  remarkable,  and  these  are  very  similar  to  the  ones  seen  in  Figures 
7  and  9. 

4.2  Boundary  Effects 

When  the  small  bead  assemblies  with  the  stress  controlled  lateral  boundary 
are  compressed  in  the  radial  direction  they  seem  to  develop  stronger  force  chains  in 
the  radial  direction  within  a  relatively  narrow  region  at  the  outer  boundary  (refer 
to  Figure  11). 

Thus,  the  small  bead  specimens  with  both  stress  and  strain  controlled  lateral 
boundaries  show  similar  behavior  when  they  have  a  comparable  fabric  (Figures  4 
and  6).  On  the  other  hand,  the  large  bead  assemblies  with  the  stress  controlled 
lateral  boundeiry  showed  larger  strain  in  the  radial  direction  as  compared  to  the 
strain  controlled  lateral  boundary  (compare  Figure  5(a)  and  (b)).  The  specimens 
with  the  two  different  lateral  boundary  conditions  are  supposed  to  have  had  a  similar 
initial  fabric.  Thus,  any  changes  in  the  isotropic  compression  test  should  be  due 
to  the  difference  in  the  boundzuy  conditions.  The  more  uniform  stress  distribution 
in  the  radial  direction  due  to  its  stress  controlled  lateral  boundary  is  thought  to 
be  responsible  for  the  larger  radial  strain  in  isotropic  compression  (see  Figure  5 

(a) )  and  large  initial  contraction  during  the  extension  tests  (see  Figures  7  (a)  and 

(b) ).  Uniformly  distributed  stress  and  the  absence  force  chains  may  cause  unstable 
internal  structures,  local  mismatches,  and  collapse,  and  may  lead  to  a  radial  larger 
strain  than  the  vertical  strain.  In  a  large  bead  assembly  there  are  not  enough 
particle  diameters  inside  the  specimen  to  develop  full  force  chains,  in  contreist  to 
the  small  bead  eissembly. 

In  the  cases  of  the  triaxial  compression  tests,  the  average  volumetric  behaviors 
of  two  specimen  with  different  particle  sizes  are  similar  to  each  other  (compare 
Figures  6(a)  and  7(a)),  both  showing  and  initial  small  densification  and  a  rather 
high  rate  of  dilation  at  a  large  strain.  The  fact  that  the  specimens  in  these  cases 
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Figure  6;  Volumetric  strains  of  small  bead  assemblies  in  triaxial  compression  and 
extension:  (a)  lateral  stress  boundary,  (b)  lateral  strain  boundary. 
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Figure  7:  Volumetric  strains  of  large  bead  assemblies  in  triaxial  compression  and 
extension  (a)  lateral  stress  boimdary,  (b)  lateral  strain  boundary. 
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(a)  Solid  large  bead  aaseably  (b)  Solid  aaall  bead  aaaeably  (c)  Hollow  saall  bead  aaaeably 


Figure  S;  Solid  and  hollow  cylinder  specimens  with  specimen  preparation  molds. 
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Figure  9:  Volumetric  strain  for  small  bead  specimens  from  hollow  cylinder  triaxial 
compression  and  extension  tests. 


(•)  Solid  lorgo  bead  asseablr  (b)  Solid  aaall  baad  aaaaablr  (c)  Hollow  saall  bead  asseablr 


Figure  11:  Schematic  diagrams  of  force  chains  developed  in  various  specimens  during 
triaxial  extension  tests. 
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had  rigid  top  and  bottom  platens  suggests  that  similar  force  chains  have  developed 
to  carry  the  major  principal  stress  across  the  sp>ecimens  in  the  vertical  direction 
(see  Figure  12). 
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Figure  12;  Schematic  diagrams  of  force  chains  developed  in  various  specimens  during 
triaxial  compression  tests. 

It  is  worth  noting  that  there  is  rather  large  scatter  in  the  data  from  the  large 
bead  assemblies  (Figure  7(a)).  This  could  be  due  to  the  insufficient  number  of 
grains  in  the  radial  directions  which  might  lead  to  a  greater  deviation  in  internal 
structure.  Another  plausible  explanation  of  the  scatter  is  the  variation  in  surface 
packing,  which  might  require  a  unique  membrane  penetration  correction  for  each 
test  case. 


5  Conclusions 

From  those  experiments,  the  researchers’  two  hypotheses  were  found  to  be  valid 
and  the  following  conclusions  can  be  made: 

1.  The  preparation  mold  creates  an  ordered  arrangement  of  particles  and  hence 
a  strong  but  localized  fabric  around  the  lateral  boundary.  The  thickness  of 
material  influenced  by  the  mold  increases  as  the  grain  size  increases.  Thus, 
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consideration  of  the  specimen  thickness  relative  to  the  grain  size  is  important; 

2.  The  strain  controlled  boundary  and  the  stress  controlled  boundary  produce 
different  volumetric  strains,  particularly  at  the  initial  stage  of  the  shearing. 
The  strain  controlled  boundary  applies  different  magnitude  of  forces  to  the 
particles  in  contact,  and  thus  generate  force  chains  through  the  specimen.  The 
formation  of  strong  force  chains  makes  the  specimen  more  expansive  from  the 
beginning.  In  contrast,  the  stress  controlled  boundary  exerts  a  uniform  forces 
on  the  particles  causing  local  compaction  and  initial  contraction  of  the  speci¬ 
men,  and 

3.  The  volumetric  behavior  of  the  thin-wall  hollow  cylinder  specimens  showed 
very  similar  behavior  to  the  solid  cylinder  specimens  with  induced  fabric  by 
preshear.  From  this  it  can  be  concluded  that  the  strong  anisotropy  in  the 
hollow  cylinder  specimen  stems  mostly  from  a  strong  initial  fabric  created  by 
the  specimen  preparation  mold. 
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Chapter  Two 


Volume  Change  in 
Inhomogeneous  Samples 


Peter  A.  Cundall 


Abstract 

Using  the  numerical  program  FLAG,  we  consider  biaxial  deforma¬ 
tions  of  an  isotropic  linearly  elastic  material  carried  out  at  fixed  means 
stress  between  combinations  of  rigid  and  flexible  boundaries.  We  show 
that  when  the  bulk  and  shear  moduli  vary  in  space  in  a  quasi-random 
fashion,  a  change  in  the  volume  of  the  sample  inevitably  occurs. 


1  Discussion 

Physical  triaxial  tests  exhibit  behavior  that  is  not  observed  in  numerical  sim¬ 
ulations  with  TRUBAL  or  in  the  results  of  simple  theoretical  calculations:  under 
constant  mean  stress,  samples  exhibit  volume  expansion  during  a  compression  test 
and  volume  contraction  during  an  extension  test.  An  explanation  is  proposed  in 
Chapter  One:  in  the  compression  test,  the  major  principal  stress  is  transmitted  via 
rigid  platens,  whereas  it  is  transmitted  through  a  flexible  membrane  in  the  exten¬ 
sion  test.  Force-chains  are  known  to  break  down  at  a  flexible  boundary,  but  are 
preserved  at  a  rigid  boundary.  Consequently,  more  “collapse” — and  thereby  volume 
decrease — occurs  in  the  extension  test. 

Although  the  proposed  explanation  is  expressed  in  terms  of  discrete  force 
chains,  a  similar  mechanism  may  occur  if  any  type  of  inhomogeneity  exists  in  a 
sample.  To  test  this,  the  program  FLAC  [Cundall  and  Board,  1988;  Cundall,  1989] 


is  used  to  simulate  biaxial  tests  on  inhomogeneous  samples,  with  combinations  of 
rigid  and  flexible  boundaries.  FLAG  models  a  continuum,  but  each  element  may  be 
assigned  a  different  set  of  properties.  In  the  tests  described  here,  an  istropic,  linear 
elastic  material  is  used — the  bulk  and  shear  moduli  vaxy  in  space  in  a  quasi-random 
fashion.  In  order  to  avoid  mesh-dependence,  the  distribution  of  moduli  is  generated 
by  a  sum  of  spacial  sine  waves:  the  shortest  wavelength  is  greater  than  the  element 
size,  and  the  longest  wavelength  is  less  than  the  sample  size  (which,  in  all  runs,  is 
10  units).  The  shear  modulus  at  any  location  (x,y)  is  computed  as 


where  (x° ,  t/°)  is  a  random  pair  of  numbers,  chosen  from  a  uniform  random  number 
distribution  in  the  range  0  <  x°  <  10  and  0  <  y°  <  10.  The  wavelength  is 
also  chosen  as  a  random  number  in  the  range  Amin  <  An  <  Amax-  For  all  tests 
reported  here,  Amin  =  2  units  and  Amax  =  10  units.  The  bulk  modulus  was  taken 
to  be  5(j/3,  corresponding  to  a  Poisson’s  ratio  of  0.25.  Each  element  in  the  grid  is 
assigned  K  and  G  according  to  its  centroid  location  (x,  y),  using  the  given  formula. 
The  constants  A  and  Go  in  equation  1  are  chosen  so  that  there  is  no  element  in 
which  G  becomes  zero  or  negative. 

Figure  1  illustrates  a  typical  distribution  of  shear  moduli,  with  magnitude 
denoted  by  intensity  of  shading;  note  that  the  biaxial  sample  is  square.  A  histogram 
of  modulus  values  for  a  similar  50x50  sample  is  provided  in  Figure  2.  Fom:  tests 
are  done  for  each  sample,  as  shown  in  Figure  3. 

Tests  (a)  and  (b)  have  rigid  platens  at  top  and  bottom;  constant  stress  is 
applied  to  the  side  boundaries.  The  tests — extension  and  compression — are  done  by 
moving  the  rigid  platens  until  the  average  platen  stresses  are  equal  and  opposite  to 
the  applied  side  stresses.  For  a  uniform  sample,  no  net  volume  change  is  observed. 
However,  in  a  non-uniform  sample  there  is  a  volume  decrease  for  the  extension 
test  and  a  volume  increase  for  the  compression  test.  For  the  case  of  an  elastic 
material,  these  two  volume  changes  are  exeictly  equal  and  opposite,  so  only  one  test 
needs  to  be  done,  for  each  boundary  condition.  In  order  to  eliminate  the  effects 
of  anisotropy  that  may  exist  in  the  sample,  two  further  tests  are  done,  with  the 
boundary  conditions  interchanged — these  correspond  to  cases  (c)  and  (d)  of  Figure 
3.  The  results  of  all  four  cases  Eire  averaged,  to  give  the  final  estimate  of  the 
volume  change.  In  order  to  achieve  the  desired  average  platen  stresses,  a  numerical 
servo-control  is  employed;  this,  and  other  sj>ecial  functions,  axe  written  in  FLAC’s 
embedded  language,  FISH.  A  sample  data  file  is  provided  in  the  Appendix. 
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Figure  1:  Shear  modulus  distribution  -  30x30  sample 


The  total  volume  change  is  found  from 


AV=  5; 


all  zones 


(<7ii  -f  (722  +  (733)  AxAy 
3K 


(2) 


where  Ax  and  Aj/  are  the  element  widths  in  x  and  y,  respectively  (the  elements 
are  rectangular).  This  formula  was  checked  against  a  direct  calculation  involving 
displacements  rather  than  stresses:  the  results  were  identical. 

In  all  runs,  a  volume  decrease  was  observed  in  the  extension  test,  and  an  equal 
volume  increase  in  the  compression  test.  Results  are  expressed  as  the  ratio,  i?,  of 
volumetric  strain  to  axial  strain.  The  following  table  records  the  values  of  R  for  the 
same  pattern  of  inhomogeneous  moduli,  but  different  grid  sizes. 


Grid 

R 

20  X  20 

2.60  X  10-2 

30  X  30 

2.74  X  10-2 

40  X  40 

2.83  X  10-2 

50  X  50 

2.90  X  10-2 
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Shear  modulus 

Figure  2:  Histogram  of  shear  modulus  values  -  50x50  sample 

It  is  seen  that  the  volume  strain  is  rather  small  -  around  3%  of  the  axial  strain. 
However,  it  is  remarkable  that  volume  strain  occurs  at  all  in  a  linear,  elastic  sample 
in  which  moduli  vary  in  a  smooth  fashion.  We  might  expect  more  of  a  “collapse” 
type  of  mechanism  in  a  discrete,  nonlinear  material,  where  extreme  local  fluctuations 
in  stiffness  axe  present  (in  the  form  of  chains  of  particles).  As  expected,  the  results 
tend  towards  an  asymptotic  value  as  the  grid  resolution  is  increased,  suggesting  that 
the  measured  volume  change  is  a  physical  effect,  rather  than  a  numerical  artifact. 

Although  the  numerical  experiments  were  intended  to  represent  (in  a  smoothed- 
out  way)  the  microscopic  fluctuations  in  moduli  of  a  granular  material,  it  is  possible 
to  view  the  inhomogeneities  as  representing  real  nonuniformities  introduced  by  the 
process  of  sample  preparation.  In  this  case,  we  might  predict  that  physical  samples 
would  exhibit  large  variations  in  volume-change  behavior,  imless  great  care  were 
taken  in  the  preparation. 
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APPENDIX;  listing  of  FLAG  data  file  for  volume-change  test 

set  log  on 
conf  extra=5 

g  30  30  ; Change  this  line,  for  different  grid  sizes 

set  echo  off 

gen  0,0  0,10  10,10  10,0 

m  e 

prop  dens=1000 

; -  Sum  of  50  special  waves;  wavelengths  picked  at  random 

; - in  the  range  L_min  to  L_max;  positions  aure  also  random 

def  cover 

g-zero  =  1 . 5e8 
top  =  jgp 
right  =  igp 
loop  i  (l,izones) 
loop  j  (l,j zones) 

shear  .mod (i,j)  =  g_zero 
end_loop 
end-loop 
loop  n  (1,50) 

xO  =  urand  ♦  10.0 
yO  =  urand  *  10.0 

lambda  =  L_min  +  urand  *  (L-mauc  -  L_min) 
fac  =  2.0  *  pi  /  lambda 
loop  i  (l,izones) 
loop  j  (l,j zones) 

XT  =  (x(i, j)+x(i+l,j)+x(i,j+l)+x(i+l,j+l))/4.0 
yr  =  (y(i,j)+y(i+l,j)+y(i,j+l)+y(i+l,j+l))/4.0 
z  =  sqrt((x0-xr)“2  +  (y0-yr)‘'2) 

shear -mod  (i,j)  =  shear-mod(i,  j)  +  amp  ♦  sinCfac  ♦  z) 
bulk_mod(i, j)  =  shear-mod(i, j)  *  5.0  /  3.0 
end-loop 
end-loop 
end-loop 

; -  check  on  lower  limit  of  G  - 

gmin  =  g-zero 
loop  i  (l,izones) 
loop  j  (l,j zones) 
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gmin  =  minCgmin.shearjnodCi,  j)) 
end^loop 
end-loop 
command 

;  -  smallest  shear  modulus  is  ... 

print  gmin 

; - should  reject  riin  if  gmin  is  negative  or  zero  ! ! 

end-Command 

end 

; -  servo  control  to  get  mean  stress  to  zero  - 

def  servo 

while-stepping 

Slim  =  0.0 

nz2  =  izones  ♦  j zones  *  2.0 

loop  i  (1, izones) 
loop  j  (l.jzones) 

sum  =  sum  +  sxx(i,j)  +  syy(i,j) 
end-loop 
end-loop 

sum  =  sum  /  nz2 
s-vel  =  min (gain*sum, 0,0001) 
if  yplat  =  1  then 
loop  i  (l,igp) 

yvel(i,l)  =  s-vel 
yvel(i,jgp)  =  -s-vel 
end-loop 
else 

loop  j  (l.jgp) 

xvel(l,j)  =  s-vel 
xvel(igp,j)  =  -s-vel 
end-loop 
end-if 
end 

; -  analysis  function  - 

def  analyse 
d-vol  =  0.0 
loop  i  (1, izones) 
loop  j  (l,j zones) 
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d_vol  =  d_vol  +  (sxxCi,  j)+syy(i,  j)+sz2(i,  j))/bulkjnod(i,  j) 
eiid_loop 
end-loop 

d_vol  =  d_vol  ♦  x(2,l)  *  y(l,2)  /  3.0 
if  yplat  =  1  then 

ax_strain  =  (ydispCl, jgp)-ydisp(l,l))  /  CyCl , jgp)-y (1 , 1)) 
else 

ax-strain  =  (xdisp(igp.l)-xdisp(l,l))  /  (x(igp,l)-x(l,l)) 
end-if 

v-strain  =  d_vol  /  (xCigp.l)  ♦  yCl.jgp)) 
rat  =  v-strain  /  ax-strain 

end 

set  echo  on 

;  Fourier  parameters  - 

set  L-min=2  L_max=10  amp=1.0e7 


cover 

set  gain=le-8 
save 

;  y-platen  tests 
set  yplat=l 
fix  y  j=l 
fix  y  j=top 
:  extension  test 
tit 

30*30  y=platen  extension 

9 

appl  p=le5  i*l 
appl  p=le5  i=right 
his  yvel  i=l  j=l 
his  sum 
step  3000 
analyse 

print  d_vol,  v-strain,  ax_strain,  rat 
save  yext.sav 

9 

rest 

:  x-platen  tests 
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set  yplat=0 
fix  X  i=l 
fix  X  i=right 
:  extension  test 
tit 

30*30  x-platen  extension 

appl  p=le5  j=l 

appl  p=le5  j=top 

his  xvel  i=l  j=l 

his  sum 

step  3000 

analyse 

print  d_vol,  v_strain,  ax_strain,  rat 

save  xext . sav 

ret 
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Chapter  Three 


Anisotropic  Elastic  Constants  of  a  Granular 

Assembly  From  Wave 
Velocity  Measurements 

Isao  Ishibashi  Tarun  K.  Agarwal 


Abstract 


Isotropic  compression,  triaxial  compression,  and  triaxial  exten¬ 
sion  tests  on  cubical  granular  specimens  (4"  x  4"  x  4")  are  performed 
in  the  laboratory.  Concurrently,  P-wave  and  S-wave  velocities  are  mea¬ 
sured  along  several  directions  relative  to  the  principal  stress  axis  in 
the  specimen.  Anisotropic  elastic  constants  are  recovered  from  those 
measurements.  The  recovered  constants  clearly  showed  the  level  of 
anisotropy  at  the  different  stages  of  the  stress  paths.  Therefore  the 
directional  wave  velocity  measurements  could  be  used  as  a  useful  tool 
to  quantitatively  identify  the  anisotropy  of  the  granular  assembly. 


1  Introduction 

The  elastic  constants  of  a  granular  assembly  and  the  wave  speeds  through  the 
media  are  fundamental  macro-properties.  These  parameters  can  also  be  used  effec¬ 
tively  for  comparison  and  vaUdation  of  numerical  and  theoretical  micromechanics 
models  through  experiments  (Jenkins  et  al.,  1988).  A  method  for  recovery  of  elastic 
constants  from  wave  speed  measurements  is  discussed  in  this  paper.  Theory  for 
wave  propagation  in  isotropic,  homogeneous  elastic  material  is  readily  available, 
however,  wave  propagation  theory  for  anisotropic  media  has  been  dealt  with  very 
little.  An  effort  is  made  to  examine  the  existing  theory,  and  its  application  to  wave 


velocity  data  obtained  from  granular  media. 


2  Elastic  Theory  of  Wave  Propagation 

Elasticity  and  inertia  are  the  properties  of  a  medium  essential  for  the  trans¬ 
mission  of  waves.  They  are  described  by  the  following  equations  in  tensor  notation: 


(a)  Equilibrium  Equation 

+  Pfi  =  P^i 

(1) 

(b)  Constitutive  Equation 

^ij  ~  Cijkl^kl 

(2) 

(c)  Strain-Displacement  Equation 

(3) 

where  (Tij  and  £ij  are  stress  and  strain  tensors,  repectively,  Cijki  is  the  fourth  order 
constitutive  tensor,  Ui  are  diplacements,  and  fi  axe  body  forces.  One  solution  of 
the  above  equations  is  the  plane  wave: 

Uk  =  Adk  exp(ia;(a:pgp  -  t)]  (4) 

where  w  is  the  real-valued  angular  frequency,  i  is  unit  imaginary  number,  qp  is 
the  component  of  the  slowness  vector  defined  as  rip/c,  with  Up  a  unit  vector  along 
the  direction  of  propagation  and  c  the  velocity  of  propagation,  Xp  are  the  coordi¬ 
nates,  and  dk  is  the  direction  of  particle  motion.  Substituting  for  the  strains  and 
accelerations  in  terms  of  the  displacement  of  a  plane  wave,  the  following  solution 
emerges: 


{CijkiqjQi  -  pSik)dk  =  0  (5) 

det  \Cijkiqjqi  -  p6ik\  =  0  (6) 

where  6ik  is  the  Kronedcor  delta  and  the  propagation  tensor  or  Green-Christoffel 
stiflFness  matrix  is  defined  by  Fife  =  CikjifijTii.  Equation  (6)  involves  a  determinant 
of  a  3  X  3  matrix,  and  defines  the  eigenvalues  that  are  related  to  three  distinct  wave 
speeds.  These  are  the  wave  speeds  of  bulk  waves  propagating  through  the  continuous 
media.  The  three  waves  are  the  pressure,  or  P,  wave  and  two  shear,  or  S  waves. 
All  three  have  mutual  orthogonal  displacement  vectors  that  are  the  eigenvectors 
of  the  problem.  In  generzd,  given  Cijki  and  the  direction  of  propagation  Uj,  the 
three  wave  velocities  and  the  associated  directions  can  be  evaluated  by  solving  the 
eigenvalue-eigenvector  problem. 
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3  Recovery  of  Elastic  Constants  from  Wave  Velocity 


The  inverse  problem  of  obtaining  elastic  constants  from  experimental  data  of 
wave  propagation  requires:  (a)  expanding  the  Equation  (6)  in  terms  of  cissociated 
elastic  constants  E.j,  which  are  related  to  the  elements  of  constitutive  tensor  Cijke', 
(b)  obtaining  wave  speeds  in  multiple  directions  in  experiments;  (c)  combining  steps 
(a)  and  (b)  for  each  wave  speed  and  developing  the  nonlinear  equation  involving 
the  direction  cosines,  wave  speeds,  elastic  constants  and  material  density;  and  (d) 
solution  of  the  multiple  nonlinear  equations  for  the  optimal  evaluation  of  the  elastic 
constants.  The  use  of  the  Newton-Raphson  procedure  with  a  good  initial  guess  gives 
an  optimal  solution.  As  described  in  Chapter  Four,  Castagnede  and  Jenkins(1989), 
and  Castagnede  et  al.  (1990)  have  applied  this  method  to  composites  and  crystals. 

The  number  of  elastic  constants  depends  upon  the  type  of  symmetry  in  the 
material  and  can  range  from  two  for  an  isotropic  material  to  21  for  a  general  elastic 
material  without  any  symmetry.  The  number  of  wave  speed  measurements  should  be 
greater  than  the  number  of  distinct  elastic  constants.  For  example,  a  transversely 
isotropic  material  described  using  five  distinct  elastic  constants,  will  require  five 
or  more  wave  speed  measurements.  Greater  redundancy  in  data  provides  higher 
acciuracy  in  the  optimized  elastic  constants. 

The  above  procedure  is  a  general  approach,  as  wave  paths  can  be  selected 
in  any  plane  in  the  continuum.  A  more  simplified  procedure  requires  that  the 
measurements  of  the  wave  vSpeeds  be  done  in  the  principal  material  planes.  For 
those  wave  paths,  the  cubic  equation  factorizes  easily  avoiding  the  complication  of 
solving  nonhnear  equations.  The  three  eigenvalues  related  to  wave  speeds  are  each 
obtained  in  sepcirate  relations  involving  the  elastic  constants,  direction  cosines  of 
the  wave  paths,  and  the  material  density.  The  difference  between  the  experimental 
velocity  and  the  velocity  from  the  elastic  constants  provides  the  error  function.  The 
summation  of  error  functions  from  veurious  wave  speed  directions  is  then  minimized 
to  obtain  an  optimal  solution  for  the  elastic  constants.  In  the  following  we  will 
discuss  the  theory  for  transversely  isotropic  materials. 

A  transversely  isotropic  stiffness  matrix  contains  5  distinct  elastic  constants. 
The  stress-strain  relation  for  such  material  is  given  as: 


'  <Tll ' 

[£^11 

Ex2 

Ei2 

0 

0 

0 

• 

'  £11  ' 

022 

Ei2 

E22 

E23 

0 

0 

0 

0 

£22 

033 

Ol2 

►  = 

E\2 

0 

E23 

0 

E22 

0 

0 

E44 

0 

0 

0 

0 

0 

£33 

2£12 

Ol3 

0 

0 

0 

E44 

0 

2£i3 

^  023  > 

.  0 

0 

0 

0 

0 

Egg 

. 

,  2£23  ' 

27 


where  Eqq  =  ^  {E22  ~  E23)  is  a  dependent  elastic  constant.  Using  the  definition  of 
the  wave  propagation  tensor  and  the  wave  path  in  the  1-2  plane,  the  components 
of  propagation  tensor  can  be  written  as; 

I'll  =  ^^i-£^ii  +nlE44;  r22  =  nlE44  +  712^^22; 

^33  = 'niE44  + -n2(E22  -  E23)  (8) 

r’12  =  (j^i2  +  £^44)711^2;  r23  =  0;  Fai  =  0  (9) 

When  equation  (6)  is  developed  using  the  above  stiffness  matrix  constants,  one 
eigenvalue  is  trivial  and  relates  to  shear  wave  velocity; 

Shear  wave;  Faa  =  pc^  =  ^J[nlE44  +  ^n^(£^22  -  ■E'23)]/p  (8) 

The  other  two  eigenvalues  are  obtained  fi-om  the  following  equation; 

(pc2)2  -  (Fu  +  F22)pc2  (FnF22  -  =  0  (9) 

The  solution  of  above  equation  gives  the  velocity  of  the  two  other  waves  as; 

P  Wave;  pcj  =  ^(-61  +  y/bl-  462);  g  ^^ve;  pch  =  ^(-61  -  -  462)  (10) 

where  61  and  62  are  given  by  the  following; 

-61  =  (Fii  -I-  F22)  =  'n\{En  -f-  044)  -I-  n2(£^22  +  -£^44);  (11) 

62  =  (F11F22  -  Ff2)  =  {n\En  +  n2£'44)  -f  (nf £^44722 £”22) 

—  (£12  +  £44)^712722  (12) 

In  transversely  isotropic  materials,  the  1-3  plane  is  similar  to  the  1-2  plane, 
hence  the  wave  velocity — elastic  constant  relations  are  similar.  Using  the  above  re¬ 
lations,  the  least  square  error  function  can  be  written  as  the  square  of  the  difference 
of  the  experimental  velocity  and  the  theoretical  velocity  given  by  the  above  equa¬ 
tions.  The  summation  of  the  errors  is  done  for  all  wave  paths.  This  least  square 
nonlinear  problem  is  then  solved  using  a  modified  Levenberg-Marqueirdt  algorithm 
and  a  finite  difference  Jacobian  firom  the  IMSL  library.  The  algorithm  requires  a 
good  initial  guess  for  the  elastic  constants,  which  is  based  on  the  expected  ratio  of 
the  various  elastic  constants. 
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4  Wave  Measurements  Experiments 

4.1  Cubic  Shear  Box 

The  generation  and  measurements  of  ultrasonic  waves  is  done  in  a  triaxial 
cubical  box  device  with  bender  bimorphs  as  transducers  (Agarwal  and  Ishibashi, 
1991)  The  specimen  is  a  4"  x  4"  x  4"  cube.  Three  independent  normal  stresses  are 
applied  through  rubber  bags  installed  inside  three  adjoining  aluminum  boundaries 
of  the  box  by  means  of  regulated  water  pressures.  The  other  three  boundaries  of 
the  box  are  made  of  solid  Incite  plates.  The  average  axial  strains  are  measured  by 
the  water  replaced  in  the  pressurized  rubber  bags.  A  glass  sphere  mixture  (two 
sizes  with  average  diameters  of  .215  mm  and  .256  mm  and  a  weight  ratio  of  1.688 
to  1)  is  used  in  the  experiments.  The  dry  mixture  is  loosely  placed  in  the  box  by 
a  spoon  and  compacted  with  a  vibrating  rod  to  the  initial  density  (approximately 
60%  of  relative  density,  equivalent  porosity  =  0.367). 

4.2  Wave  Path  Configurations 

All  wave  paths  are  located  at  the  central  section  of  the  cube  on  the  vertical 
planes  as  shown  in  Figure  1. 


Acoustic  Transducers 


Figure  1:  Wave  Path  Configuration 

Overall,  six  wave  paths  for  compression  and  three  wave  paths  for  shear  are 
employed.  The  number  of  wave  paths  wjis  selected  based  upon  the  five  dynamic 
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elastic  constants  required  for  transversely  isotropic  granular  assemblies.  The  shear 
waves  propagate  along  three  central  axes  of  the  cube. 

4.3  Stress  Paths 

A  specimen  of  granular  material  possess  an  initial  anisotropy  that  is  due  to 
its  deposition  under  gravity.  This  anisotropy  is  considered  here  to  be  transversely 
isotropic.  This  symmetry  is  also  assumed  for  all  applied  stress  paths  as  the  paths 
are  irrotational.  In  Test  one,  an  isotropic  stress  is  increased  from  4  psi  to  20  psi 
in  increments  of  4  psi.  In  Test  two,  after  an  isotropic  consolidation  of  20  psi,  a 
compression  test  is  conducted  until  failure  maintaining  a  constant  mean  stress  (20 
psi).  The  vertical  (axis  1)  normal  stress  is  increased  while  the  normal  stresses  in 
two  horizontal  directions  (axes  2  and  3)  are  decreased  equally.  In  Test  three,  after 
an  isotropic  consolidation  of  20  psi,  the  specimen  is  sheared  in  extension  (a  decrease 
in  the  vertical  stress  and  equal  increase  in  horizontal  stresses),  again  maintaining  a 
constant  mean  stress  of  20  psi. 

5  Results  and  Discussions 

Wave  speed  measurements  for  P-waves  (6  paths)  and  S-waves  (2  paths)  are 
presented  in  Tables  I,  II,  and  III.  The  angle  (5  refers  to  the  inclination  of  the  wave 
path  from  the  vertical  in  the  vertical  planes  of  propagation.  The  first  five  wave 
speed  measurements  are  conducted  in  the  1-2  plane  and  the  sixth  wave  path  lies  in 
the  1-3  plane.  Both  shear  wave  paths  he  in  1-2  plane. 

Table  I  shows  that  in  all  the  measurements  the  wave  velocities  increase  with 
increasing  confining  pressure.  Different  wave  velocity  measurements  in  the  two 
orthogonal  directions  =  0°  and  =  90°  indicate  the  existence  of  an  inherent  fabric 
formed  during  the  deposition  of  the  specimen.  The  strain  under  isotropic  pressure 
showed  that  the  specimen  was  more  compressible  in  the  horizontal  direction  than 
in  the  vertical  direction  as  seen  in  Figure  2. 

This  is  expected,  since  the  deposition  under  gravity  will  induce  a  relatively 
Ijirger  number  of  contacts  in  the  vertical  direction  and  thus  the  specimen  will  be 
less  compressible  in  that  direction.  In  terms  of  the  static  elastic  modulus,  then, 
the  specimen  was  softer  in  the  horizontal  direction  than  in  the  vertical  direction. 
Similarly,  wave  speeds  were  larger  in  the  vertical  direction  compared  to  those  in 
the  horizontal  direction.  It  is  noted,  however,  that  Knox  et  al.  (1982)  conducted 
P-wave  tests  on  the  cubic  specimens  of  sand  that  showed  a  higher  wave  velocity 
in  the  horizontal  directions.  It  is  difficult  to  explain  their  experimental  data  from 
consideration  of  the  depositional  fabric. 

Table  II  shows  wave  velocities  until  near  failure  in  the  compression  test.  Upon 
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Figure  2:  Strain  responses  during  istoropic  pressure  test. 

the  application  of  12  psi  maximum  shear  stress,  the  increase  in  the  P-wave  velocity  in 
the  vertical  direction  is  5.4%  while  the  decrease  in  P-wave  velocity  in  the  horizontal 
direction  is  between  13.4%  and  10.7%  .  The  two  shear  wave  meeisurements  show  very 
small  changes  (a  5  to  6%  drop).  This  is  significantly  different  from  the  decreases  in 
the  shear  wave  velocity  of  15%  to  20%  in  the  hollow  cylindrical  experiments  (Chen 
et  al.,  1988). 

Table  III  summarizes  wave  velocities  for  the  extension  test.  The  P-wave  velocity 
drops  by  16.9%  in  the  vertical  direction  and  increases  by  9.7%  in  the  horizontal 
direction  during  the  shear.  The  shear  wave  velocity  also  shows  greater  changes 
compared  to  the  compression  test. 

Table  IV  and  Figure  3  show  recovered  elastic  constants  for  the  transversely 
isotropic  specimen  during  the  isotropic  stress  application. 

Both  dynamic  elastic  and  shear  modulus  are  larger  in  the  1-2  vertical  plane 
than  in  the  2-3  horizontal  plane.  The  changes  in  elastic  modulus  in  two  directions 
follow  an  often- used  equation  E  =  Ka°,  where  K  and  a  are  constants  and  a  is  the 
isotropic  confining  stress.  Values  of  K  and  a  for  different  directions  are  calculated 
as  follows:  for  En:  K  =  351.0,  a  =  .169;  for  S22:  K  =  225.1,  a  =  .255;  for  E^: 
K  =  87.7,  a  =  .296;  for  E^%-.  K  =  79.6,  a  =  .325.  Tests  conducted  in  hollow 
cylindrical  samples  (Chen  et  al.,  1988)  showed  a  =  .40  for  the  shear  modulus  in  the 
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Table  I:  Wave  velocity  in  m/sec  during  isotropic  stress  application. 


Pres¬ 

sure 

Pu- 

Wave 

Pn- 

Wave 

Pn- 

Wave 

Pu- 

Wave 

Pn- 

Wave 

P.S- 

Wave 

Sl2- 

Wave 

Sn- 

Wave 

(psi) 

at 

at 

at 

at 

at 

at 

at 

at 

p-o® 

P-21® 

P-46®  i 

P-70® 

P-90® 

P-90® 

P-90® 

P-0® 

4 

432.2 

425.8 

451.1 

425.1 

415.7 

435.7 

246.4 

240.2 

8 

455o 

430.9 

465.6 

429.4 

449.7 

473.0 

268.7 

263.7 

12 

481.4 

461.5 

489.1 

487.1 

500.6 

503.4 

291.7 

279.4 

16 

491.7 

486.6 

527.6 

520.0 

524.5 

524.6 

314.4 

295.9 

20 

510.2 

500.1 

544.1 

538.1 

556.0 

538.6 

298.3 

323.8 

Table  11;  Wave  velocity  in  m/sec  during  compression  test. 


Shear 

Stress 

Pu- 

Wave 

Pu- 

Wave 

Pi2- 

Wave 

Pu- 

Wave 

Pa- 

Wave 

Pa- 

Wave 

Sl2- 

Wave 

^12- 

Wave 

(psi) 

at 

at 

at 

at 

at 

at 

at 

at 

P-0® 

P-21® 

P-46® 

P-70® 

P-90® 

P-90® 

P-90® 

P-0® 

0 

510.2 

500.1 

544.1 

538.1 

556.0 

538.6 

2983 

323.8 

1.5 

517.1 

503.5 

540.6 

528.8 

537.8 

533.1 

298.4 

325.1 

3.0 

512.1 

497.6 

539.8 

5173 

525.9 

520.8 

2923 

318.8 

4J 

5253 

5043 

524.7 

490.7 

547.0 

530.9 

274.8 

322.1 

6.0 

535.1 

498.8 

521.9 

530.6 

5303 

510.7 

2903 

324.2 

7J 

540.7 

505.1 

543.4 

497.1 

508.8 

493.6 

274.8 

323.6 

9.0 

5493 

500.6 

526.6 

480.0 

491.9 

455.9 

274.1 

3193 

10.5 

5643 

502.8 

520.1 

460.8 

477.0 

423.4 

2923 

3193 

12.0 

569.6 

5053 

507.7 

448.0 

4643 

413.4 

292.8 

319.9 

Table  HI:  Wave  velocity  in  m/sec  during  extension  test. 


Shear 

Stress 

(psi) 

Pu- 

Wave 

at 

P-0® 

Pa- 

Wave 

at 

P-21" 

Pa- 

Wave 

at 

P-46® 

Pa- 

Wave 

at 

P-70® 

Pa- 

Wave 

at 

P-90® 

Pa- 

Wave 

at 

P-90® 

Su- 

Wave 

at 

p-90® 

Su- 

Wave 

at 

P-0® 

90 

517.4 

4953 

5173 

500.8 

534.4 

538.6 

2973 

323.6 

13 

489.1 

4803 

546.1 

4923 

5323 

533.1 

3023 

314.9 

— 

471.8 

4653 

5123 

541.7 

550.9 

5539 

306.4 

310.7 

43 

4653 

4483 

522.7 

565.7 

617.0 

619.1 

311.8 

311.1 

mm 

4393 

422.8 

509.9 

575.6 

6373 

634.7 

3123 

300.1 

\  73 

410.7 

4003 

499.4  ' 

585.1 

652.8 

6536 

308.8 

290.6 

8.6 

3833 

366.6 

496.6 

533.1 

652.9 

653.1 

309.1 

2933 
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Table  IV:  Dynamic  elastic  constants  in  MPa  dtiiing  isotropic  stress  application. 


Pressure 

(psi) 

Ell 

Ei2 

E- 

E« 

E« 

4 

286.6 

149.5 

282.4 

95J 

90.6 

8 

313.7 

93.9 

113.4 

109.2 

12 

352J 

104.1 

393.0 

133.6 

122.6 

16 

370.1 

166.8 

430.8 

155.2 

137.4 

20 

397.0 

129.0 

467.7 

139.7 

164.6 

Table  V:  Dynamic  elastic  constants  in  MPa  during  compression  test. 


Shear 

Stress 

(psi) 

Eu 

Ei2 

E-2 

E« 

Efi, 

0 

397.0 

129.0 

467.7 

139.7 

164.6 

1.5 

407.4 

1223 

448.4 

139.8 

165.8 

3.0 

398.0 

140.2 

4273 

134.1 

1593 

4_5 

425.1 

55.9 

438.6 

1183 

162.8 

6.0 

427.0 

77.8 

4323 

132.5 

165.0 

7.5 

435.9 

1273 

390.6 

1183 

164.4 

9.0 

444.2 

103.9 

351.6 

118.0 

160.1 

103 

463.8 

83.2 

3173 

1343 

1603 

12.0 

473.3 

49.7 

301.4 

134.6 

160.7 

Table  VI:  Dynamic  elastic  constants  in  MPa  during  e.ttension  test. 


Shear 

Stress 

(psi) 

E,i 

Eu 

E« 

E« 

0 

4103 

45.1 

440.7 

138.8 

164.4 

13 

364.6 

140.4 

4303 

143.6 

155.6 

341.0 

763 

478.9 

1473 

151.6 

43 

330.7 

0.0 

585.1 

152.6 

151.9 

6.0 

2933 

0.0 

6163 

1533 

135.4 

73 

256.7 

0.0 

639.7 

149.7 

123.0 

8.6 

2173 

0.0 

638.1 

150.0 

118.7 

Figure  3:  Dynamic  elastic  constants  during  isotropic  pressure, 
horizontal  plane. 

The  elastic  modulus  E12  is  found  to  change  erratically  with  change  in  the  pres¬ 
sure.  The  other  elastic  moduli  follow  much  more  systematic  chages  in  magnitude. 
The  exact  identification  of  the  shear  wave  arrivaJ  was  difficult  in  some  cases  and 
this  might  have  created  some  errors  in  Eu. 

Table  V  shows  elastic  modulus  during  the  compression  test  and  Figure  4  plots 
the  same  data  with  respect  to  the  maximum  shear  stress. 

The  increase  in  En  is  smaller  than  the  decrease  in  £^22-  This  may  indicate  that 
the  increase  in  number  of  contact  normals  in  the  vertical  direction  is  low  compared 
to  the  loss  in  number  of  contact  normals  in  and  eiround  the  horizontal  plane.  Note 
that  overall  contacts  are  lost  during  shearing  in  numerical  tests  (Ishibashi  et  al., 
1989) .  These  fabric  changes  are  reflected  only  by  moderate  gains  in  elastic  modulus 
in  vertical  direction  (10.6%)  and  significant  drops  in  horizontal  modulus  (29%). 
Therefore  the  dynamic  elastic  modulus  is  an  important  reflection  of  fabric  changes 
inside  the  granular  material  diuing  shear.  The  changes  in  shear  modulus  axe  10% 
to  17%.  In  comparison,  the  tests  with  hollow  cylindrical  specimens  has  shown 
significant  drops  in  shear  modulus  during  shear.  The  difference  in  the  wavelengths 
in  the  two  type  of  tests  could  be  a  reason  for  this  difference.  In  the  cube  specimens, 
the  wavelength  is  about  10  particle  diameters  while  the  hollow  cylindrical  tests  use 
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ts.  MPa 


Figlire  4:  Dynamic  elastic  constants  during  compression  test. 


Figure  5:  Dynamic  elastic  constants  during  extension  test. 
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a  wavelength  of  the  order  of  specimen  height  (7.6  inches). 

Table  VI  and  Figure  5  show  changes  in  elastic  modulus  during  the  extension 
tests.  The  changes  in  the  two  directional  elastic  modulus  E\\  and  E22  are  relatively 
large  compared  to  the  compression  test.  The  specimen,  initially  has  less  contact 
normals  in  the  horizontal  direction  compared  to  the  vertical  direction  due  to  depo¬ 
sition  under  gravity.  During  the  extension  test,  increases  in  the  contact  normals  in 
the  horizontal  direction  axe  significantly  larger  than  the  increases  of  contact  nor¬ 
mals  in  the  vertical  direction  during  the  compression  tests.  This  difference  probably 
causes  large  increase  in  £^22-  Simultaneously,  the  contact  normals  decrease  in  the 
vertical  direction  and  thus  the  ela.stic  modulus  E\\  reduces.  The  shear  modulus  in 
the  horizontal  plane  shows  reduction  similar  to  the  compression  test.  It  should  be 
noted  that  the  initial  isotropic  stage  of  two  specimens  (Figures  4  and  5)  are  different 
and  thus  caused  slightly  dissimilar  En  and  E22  values. 

Conclusion 

A  general  procedure  for  relating  dynamic  elastic  constants  of  anisotropic  mate¬ 
rials  to  wave  speeds  is  presented  and  detailed  mathematical  relations  for  transvesely 
isotropic  specimens  has  been  discussed.  The  approach  is  applied  successfully  to  the 
multi  directional  wave  speed  measurements  in  cubic  specimens  of  granular  assembly 
of  spherical  particles  for  three  different  stress  paths.  In  an  isotropic  stress  test,  the 
compression  wave  speeds  in  horizontal  direction  were  lower  compared  to  the  speeds 
in  the  vertical  direction  due  to  depositional  fabric.  The  recovered  dynamic  elas¬ 
tic  moduli  also  followed  this  pattern.  This  behavior  is  similar  to  the  static  elastic 
modulus  as  is  calculated  from  the  three  axial  strain  measurements.  The  recovered 
elastic  moduli  showed  a  systematic  change  during  compression  and  extension  tests. 
The  importance  of  changes  in  internal  structure  as  characterized  by  contact  normals 
is  quahtatively  associated  with  the  recovered  elastic  moduli.  In  sununary,  the  eval¬ 
uation  of  dynamic  elastic  moduli  for  granular  assembly  is  a  useful  tool  for  studying 
internal  character  of  granular  materials. 
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Chapter  Four 


Determination  of  the  Principal  Acoustic  Axes 

in  Anisotropic  Solids 
from  Wavespeed  Measurements 

B.  Castagnede  J.T.  Jenkins 


Abstract 

We  outline  and  provide  references  to  procedures  that  are  em¬ 
ployed  to  determine  the  anisotropy  of  linearly  elastic  solids  from  ex¬ 
perimental  measurements  of  the  speeds  of  elastic  waves  propagating 
through  them  and  apply  one  of  the  numerical  methods  to  the  analysis 
of  some  data.  We  indicate  how  the  procedures  might  be  extended  t  > 
characterize  the  anisotropy  of  a  granular  material. 


1  Disscussion 

Prediction  of  the  incremental  behavior  of  a  homogeneous  sample  of  a  granu¬ 
lar  material  requires  information  about  its  current  state.  The  identification  of  the 
variables  that  characterize  the  state  of  the  materials  is  currently  an  area  of  active 
research,  e.g.  (1).  Certain  to  be  included  on  any  list  of  state  variables  are  the 
porosity  and  stress.  Other  candidates  are  the  orientational  distribution  of  contacts 
and  contact  area.  The  hope  is  that  measurements  of  wavespeed  and  attenuation 
in  relatively  homogeneous  samples  can  help  to  characterize  the  state  subsequent  to 
some  loading  history.  For  example,  Chen,  et  al.  (2,3)  measure  the  speed  of  a  shear 
wave  propagating  down  the  axis  of  a  hollow  cylindrical  sample  in  a  torsional/triaxial 
device.  They  relate  the  observed  decrease  in  wavespeed  with  increasing  deviatoric 
stress  to  a  net  loss  of  contacts  associated  with  the  development  of  anisotropy.  Here 


we  wish  to  discuss  how  wavespeed  measurements  in  classical  linear  elastic  solids  can 
provide  rather  complete  characterizations  of  the  acoustic  anisotropy.  Our  expecta¬ 
tion  is  that  similar  methods  can  be  employed  to  determine  the  evolving  anisotropy 
in  granular  materials. 

An  elastic  sohd  is  typically  produced  with  a  symmetry  that  does  not  change 
when  the  material  is  subjected  to  small  strains.  An  example  is  an  epoxy  matrix 
reinforced  by  parallel  fibers.  The  linear  elastic  behavior  is  governed  by  the  stiffness 
tensor  Cijki  that  relates  the  stress  Uj  and  the  infinitesimal  strain  eki  ■ 

^ij  ~  Cijkl&kl-  (1) 

Because  of  the  symmetry  of  the  stress  and  the  strain,  the  stiffness  tensor  satisfies 
the  conditions  Cijki  =  CjM  =  Cijik-  The  existence  of  a  strain  energy  provides  the 
additional  condition  Cijki  =  CkUj-  Other  reductions  in  the  number  of  independent 
components  of  the  stiffness,  or  the  elastic  constants,  depend  upon  the  symmetry  of 
the  material  (4). 

If  the  solid  is  subjected  to  a  large  deformation,  its  symmetry  is,  in  general, 
changed.  For  example,  if  the  epoxy  composite  is  greatly  compressed  in  a  direc¬ 
tion  perpendicular  to  the  fibers,  its  synunetry  is  transformed  firom  hexagonal  to 
orthorhombic.  The  description  of  wave  propagation  through  a  material  that  is  ex¬ 
periencing  a  large  deformation  is  complicated  by  the  presence  of  the  stress  t^j  that 
maintains  the  deformation.  In  this  case  the  increment  t'^j  of  stress  associated  with 
the  wave  results  from  both  the  additional  increment  of  strain  Cij  and  the  increment 
in  the  antisymmetric  part  Wij  of  the  displacement  gradients, 

t[j  =  Wikt^j  —  t^k'^kj  +  (l/2)(eiA:<fcj  +  tik^kj)  —  ^kkt^j  +  CijklCkl,  (2) 

where  the  stiffness  tensor  Cijki  depends  upon  the  symmetry  of  the  finitely  deformed 
material  (5). 

The  situation  in  granular  materials  is  complicated  by  the  fact  that  a  part  of  the 
anisotropy  that  is  induced  by  a  history  of  deviatoric  stress  remains  upon  unloading. 
However  it  is  precisely  this  anisotropy  that  endows  the  material  with  the  memory 
of  its  history.  The  difficulty  is  that  this  anisotropy  does  not,  in  general,  fall  into  one 
of  the  common  classes  (6).  Nevertheless  we  think  that  it  is  worthwhile  to  review 
existing  schemes  for  the  determination  of  the  elastic  constants  and  the  orientation  of 
the  principal  acoustic  axes  in  crysteds  and  in  anisotropic  engineering  materials.  This 
is  a  very  practical  problem  that  has  been  st\idied  and  solved  for  various  material 
synunetries  using  numerous  numerical  methods. 

When  the  symmetry  of  the  material  is  known,  the  symmetry  axes  induce  a 
natural  orthonormal  basis.  The  directions  peirallel  to  these  vectors  are  called  the 
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principal  acoustic  ajces.  The  elastic  constants  are  to  be  determined,  so  we  need  to 
write  down  the  characteristic  equation  associated  with  propagation  in  a  direction 
with  cosines  rii, 


det  iFij  -  pv^Sijl  =  0-  (3) 

Here  F ij  is  the  Green-Christoffel  tensor,  given  by 

Fife  —  CijklTljTll,  (4) 

p  is  the  density  of  the  material,  6ij  is  the  Kronecker  delta,  and  v  is  the  wavespeed 
of  a  bulk  mode. 

When  the  acoustic  axes  of  the  material  are  not  known,  these  equations  are  still 
valid,  but  the  natiaral  basis  is  unavailable.  The  direction  cosines  n[  are  then  given 
with  respect  to  some  other  basis  convenient  to  the  experiment: 

Tlj  —  (lijTlj^  (5) 

where  atj,the  transformation  matrix,  is  given  in  terms  of  Euler  angles  (7)  that  must 
be  determined  as  part  of  the  solution  to  the  problem. 

Equation  (3)  is  a  cubic  equation  for  the  eigenvalues  of  the  propagation  tensor, 
pv^.  These  eigenvalues  are  measured  experimentally  for  a  number  of  orientations 
n'.  Then,  in  order  to  obtain  the  elastic  constants  of  the  material,  equation  (3)  must 
be  inverted. 

To  accomplish  this  inversion,  several  methods  have  been  proposed  during  the 
last  40  years  including: 

A  numerical  procedure  (8)  and  series  expansion  (9); 

Approximations  valid  near  the  principal  axes  of  symmetry  (10); 

Perturbation  series  expansion  with  an  iterative  numerical  method  (11-15); 

Use  of  the  three  invariants  of  the  ChristoflFel  matrix  (16,  17);  and 

Algebraic  manipulations  of  the  characteristic  equation  (18,  19). 

In  the  general  case,  the  use  of  this  last  method  provides  a  new  cubic  equa¬ 
tion  (20,  21)  for  the  determination  of  the  colunm  vector  Xa  of  unknown  elastic 
constants. 


{1/3\)AabcXaXbXc  +  il/2\)BABXAXB  +  CaXa  +  D=0  (6) 

where  the  dimension  of  Xa  ranges  from  3  for  cubic  symmetry  to  21  in  triclinic 
materials.  The  coefficients  depend  upon  the  density,  the  direction  cosines,  and  the 
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Euler  angles.  They  eire  given  elsewhere  (21)  for  the  hexagonal  and  orthohombic 
symmetries. 

Generally,  the  systems  of  equations  (6)  is  over  deter  mined,  as  there  is  more 
experimental  data  than  unknown  elastic  constants.  An  optimal  solution  of  the 
overdetermined  system  is  obtained  using  a  Newton- Rnphson  numerical  scheme.  The 
recovery  of  the  Euler  angles  has  been  achieved  using  this  method  or  other  procedures 
for  various  systems  of  symmetry.  The  problem  has  been  solved  for  the  cubic  (17, 
22),  the  hexagonal  (21)  and  even  for  the  more  intricate  monoclinic  (23,  24)  systems. 

When  the  type  of  symmetry  is  not  known,  equation  (6)  may  still  be  employed 
with  the  dimension  of  Xa  taken  to  be  the  number  of  independent  elastic  constants 
of  a  candidate  symmetry.  If  the  material  exhibits  more  symmetry,  the  calculated 
stiffness  tensor  will  be  seen  to  have  fewer  independent  components  than  initially 
expected  (25).  Using  such  a  procedure,  experimental  data  on  some  engineering 
materials  may  be  fit  equally  well  by  either  of  two  symmetries.  In  granular  materials 
we  expect  this  to  be  the  case  more  often  than  not,  but  we  anticipate  greater  difficulty 
in  recognizing  what  symmetry  we  may  be  close  to. 

2  An  Example 

The  above  approach  is  illustrated  in  one  of  the  simplest  cases,  the  recovery  of  a 
single  Euler  angle  for  the  hexagonal  system  of  symmetry.  This,  for  example,  is  the 
case  when  considering  the  recovery  of  the  ahgnment  of  the  fibers  in  a  unidirectional 
composite  material.  Figure  1  provides  the  geometry  for  this  problem  when  dealing 
with  a  thin  plate  of  the  composite.  We  suppose  that  this  plate  was  cut  with  an 
angular  parallax  $23  and  that  this  parameter  is  unknown.  The  idea  is  then  to 
propagate  eleistic  waves  in  the  planes  (1,2')  and  (1,3')  and  assume  them  to  be 
principal  planes.  In  fact,  the  principal  planes  are  the  coordinate  planes  of  the 
unprimed  coordinate  system. 

Figure  2  shows  the  three  dimensional  mapping  of  the  slownesses  (the  inverse 
of  the  wave  speeds)  for  the  quasi-trjmsverse  (q.T.)  mode  of  propagation  in  the 
hexagonal  symmetry.  In  general,  there  are  two  other  bulk  waves  that  should  be 
considered.  These  are  omitted  here  in  order  to  simplify  the  drawing.  This  mapping 
is  a  surface  of  revolution  that  intersects  the  (1,2)  plane  in  a  circle.  The  3-axis  is  the 
true  orientation  of  the  fibers  of  the  composite.  However  these  fibers  are  believed  to 
be  oriented  along  the  axis  3'  in  the  coordinate  system  of  the  experiment. 

When  the  elastic  waves  zire  propagated  in  the  plane  (1,3),  the  two  dimensional 
slice  (TE)  is  experimentally  generated.  Due  to  the  invariance  of  the  slownesses  to 
a  rotation  around  axis  3,  other  slices  such  as  {B\E)  and  {B2E)  are  identical  to 
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Figure  3:  Normalized  slowness  for  a  imidirectional  glass/epoxy  composite  material 
in  the  istropic  plane  (1,2).  Crosses:  Experimental  slownesses;  Continuous  lines: 
Recovered  slowness  without  angular  parallax. 


2.2' 


Figure  4:  Noramalized  slowness  for  a  unidirectionzd  glass/epoxy  composite  mate¬ 
rial  in  the  anisotropic  plane.  Crosses:  Experimental  slowness;  Continuous  lines: 
Recovered  slowness  without  angular  parrallax-plane  (1,3);  Dashed  lines:  Recovered 
slowness  with  8  degrees  paxallax-plane(l,3'). 

segment  (ET).  By  introducing  the  Euler  angle  023  >  we  determine  what  the  situation 
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is  when  the  propagation  taJces  place  in  (1, 3')  plane.  In  this  case,  we  have  to  examine 
the  projection  of  the  straight  line  {D)  onto  the  surface  of  revolution,  i.e.  the  curve 
(C).  This  new  slowness  curve  crosses  (C/)  and  (C//),  indicating  that  it  is  diflFerent 
from  the  corresponding  curve  in  the  principal  plane  (1,3). 

Figure  3  and  4  show  results  for  a  unidirectional  fiberglass/epoxy  composite.  In 
Figiure  3  there  is  isotropy  in  the  plane  (1,2)  that  is  orthogonal  to  the  fibers.  In  Fig¬ 
ure  4  the  slownesses  in  the  (1,3)  and  (1,3')  planes  are  given.  The  crosses  represent 
experimental  data  obtained  by  using  an  advanced  ultrasonic  spectro-interferometer 
(19).  The  lines  axe  the  recovered  slownesses  generated  after  optimally  solving  equa¬ 
tion  (6).  The  solid  lines  correspond  to  ^23  =  0°  (no  angular  parallax)  and  the  dashed 
hnes  are  for  the  case  023  =  8°  (the  sample  is  simply  rotated  on  its  stage  by  8°  around 
the  1-axis).  The  shape  of  the  recovered  slownesses  with  and  without  parallax  is  sim¬ 
ilar.  It  is  interesting  to  note  that  there  is  no  change  in  the  quasi-longitudinal  {q.L.) 
wavespeed  along  the  1-axis  (the  point  labelled  I  in  Figure  4)  due  to  the  fact  that 
1  =  1'.  For  the  same  mode,  the  wavespeed  along  the  3-axis  (the  point  labelled  J) 
is  smaller  when  023  =  8°  compared  to  the  case  with  no  parallax.  Again,  this  is  not 
surprising  as  the  largest  wavespeed  of  this  mode  is  along  the  fibers.  After  rotating 
the  sample  by  023,  J  no  longer  corresponds  to  the  orientation  of  the  fibers. 

Using  a  double-iterative  algorithm  derived  from  the  set  of  equations  (6)  in  the 
primed  coordinate  system  (21),  we  were  able  to  recover  a  value  for  the  Euler  angle  of 
8.368°  ±0.319°.  The  elastic  constants  may  be  recovered  with  even  greater  accuracy. 
The  precision  in  the  recovery  weis  obtained  using  a  numerical  simulation  in  which 
the  parameters  were  selected  to  closely  match  those  of  the  experiment  (21). 
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Chapter  Five 


Anisotropic  Elasticity  for  Random 
Arrays  of  Identical  Spheres 

James  T.  Jenkins 


Abstract 


We  derive  a  relation  betv/een  the  average  stress  and  the  aver¬ 
age  strain  in  a  homogeneous  deformation  of  a  transversely  isotropic 
but  otherwise  random  array  of  identical  spheres  that  interact  through 
noncentral  contact  forces.  We  assume  that  the  displacement  of  a  con¬ 
tact  relative  to  a  center  may  be  calculated  from  the  average  strain 
of  the  aggregate  and  the  average  rotation  of  the  spheres  relative  to 
that  of  the  sample.  In  a  homogeneous  deformation  of  an  anisotropic 
aggregate,  the  average  rotation  of  the  spheres  is  determined  from  the 
requirement  that  the  stress  be  symmetric.  In  general,  it  differs  from 
the  average  rotation  of  the  sample.  Then,  having  raised  the  issue  of 
the  different  rotations,  we  focus  on  strongly  inhomogeneous  deforma¬ 
tions,  derive  general  expressions  for  the  stress  and  couple  stress,  and, 
as  an  illustration,  calculate  the  constitutive  relations  for  small  strains 
and  rotations  of  an  isotropic  array  with  elastic  contacts. 


1  Introduction 

We  are  interested  first  in  the  average  small  strain  response  of  random  granu¬ 
lar  assemblies  that,  by  virtue  of  the  process  by  which  they  are  formed,  possess  an 
inherent  anisotropy  [1].  This  interest  originated  in  attempts  to  explain  the  volume 
changes  observed  in  trieixial  extension  and  trieixial  compression  of  a  granular  mate¬ 
rial  at  a  fixed  pressure.  These  experiments  were  carried  out  on  hollow  cylindrical 
samples  of  essentially  identical  glass  spheres  confined  within  a  membrane  in  a  true 


triaxial/torsion  device  [2].  As  Jenkins  [3]  shows,  the  observed  volume  changes  can 
be  explained  in  the  context  of  the  anisotropic  elasticity  theory  developed  here  for 
homogeneous  deformations.  However,  as  described  in  Chapter  1,  more  recent  ef¬ 
forts  to  simulate  these  volume  changes  numerically  have  indicated  that  boimdaries 
probably  played  the  dominant  role  in  the  experiments.  In  any  case,  the  theory  for 
homogeneous  deformations  of  an  anisotropic  aggregate  has  one  important  feature 
with  several  interesting  consequences.  Because  the  contact  forces  are  noncentral, 
symmetry  of  the  stress  requires  that  the  average  rotation  of  the  particles  be  differ¬ 
ent  from  the  average  rotation  of  the  sample.  We  illustrate  this  in  the  context  of  a 
theory  derived  for  arrays  in  which  the  geometric  arrangement  exhibits  transverse 
anisotropy  and  the  strains  are  so  small  that  the  behavior  is  elastic. 

The  interest  in  strongly  inhomogeneous  deformations  arises  from  recent  at¬ 
tempts  to  explain  the  structure  of  narrow  regions  of  rather  extreme  inhomogeneity, 
called  shear  bands,  in  granular  materials.  An  example  of  this  activity  is  the  work 
of  Miihlhaus  and  Vardoulakis  [4].  Such  attempts  involve  the  introduction  of  the  ad¬ 
ditional  degrees  of  freedom  associated  with  the  independence  of  the  mean  rotation 
of  the  particles.  In  existing  treatments,  energy  arguments  are  used  to  obtain  the 
constitutive  relations  that  link  the  stress  and  couple  stress  to  the  strain  and  rota¬ 
tions.  Here,  we  illustrate  how  such  relations  may  be  obtained  from  considerations 
of  force  and  moment  alone  and  provide  a  simple  example  that  is  appropriate  for  the 
small  strain  elastic  response  of  an  isotropic  assembly.  In  these  developments,  the 
anisotropy  of  the  aggregate  is  not  crucial,  although  in  granular  materials  it  is  likely 
that  some  form  of  anisotropy  is  always  present. 

For  simplicity,  we  focus  our  attention  on  an  idealized  material  consisting  of 
identical  spherical  grains.  We  suppose  that  the  particles  are  in  contact  and  that 
the  location  of  their  centers  is  random.  In  reality,  arrays  of  identical  spheres  often 
organize  themselves  into  domains  in  which  their  centers  are  arranged  periodically, 
but  we  ignore  this  or  imagine  that  there  is  sufficient  dispersion  in  the  diameters  to 
prevent  the  formation  of  such  crystalline  regions. 

2  Homogeneous  Deformations 

We  suppose  that  the  diameter  of  a  sphere  is  and  that  there  are  n  of  them  per 
imit  volume.  We  take  a  to  be  the  unit  vector  from  the  center  of  a  sphere  to  a 
contact  on  its  surface  and  introduce  the  orientational  distribution  of  contacts  A{a ) 
so  that  A{a  )da  is  the  probable  number  of  contacts  in  the  element  da  of  solid  angle 
centered  at  a  .  In  homogeneous  situations,  this  distribution  fimction  is  independent 
of  position  and,  because  a  contact  is  common  to  two  spheres,  A{—a)  =  A(a). 
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When  the  distribution  of  contacts  is  isotropic,  A{a)  =  k/Air,  where  k,  called  the 
coordination  number,  is  the  average  number  of  contacts  per  particle.  When  the 
contact  distribution  is  not  isotropic,  but  there  are  orthogonal  planes  of  symmetry 
or  an  axis  of  symmetry,  the  distribution  of  contacts  can  be  expressed  in  terms  of  a 
symmetric,  traceless  tensor  A  [5,  6]  : 


k 

A{ol)  = AijUiUj).  (1) 

For  transversely  isotropic  materials,  A  may  be  expressed  in  terms  of  the  unit  vector 
h  in  the  direction  of  the  axis  of  anisotropy  and  the  strength  £  (1)  of  the  anisotropy: 

Aij  =  -e{6ij  -  Zhihj),  (2) 


where  —  ^  <  £  <  1.  Then 


X(a)  =  |j:[(l-e)  +  3£(fciaifl.  (3) 

The  force  F{a )  exerted  by  the  sphere  at  a  contact  with  orientation  (3)  has 
components  parallel  and  perpendicular  to  a : 


Fi  =  Poi  -  Ti,  (4) 

where  T  a  =0.  In  a  homogeneous  deformation,  we  assume  that  Fis  independent  of 
position.  Because  the  force  exerted  by  the  sphere  at  a  contact  is  equal  and  opposite 
to  the  force  exerted  on  the  sphere,  F{—a )  =  — ).  The  normal  and  tangential 
components  depend  upon  the  displacement  of  a  contact  relative  to  the  center  and 
are,  in  general,  functions  of  (4).  Following  Hertz  [7],  the  normal  component  P  is 
related  to  the  magnitude  (4)  of  the  normal  component  of  the  displacement  s: 

P  =  (5) 


where  M  is  given  in  terms  of  the  shear  modulus  fj,  and  Poissons  ratio  v  of  the 
material  of  the  spheres  by 

9,/3(l-i/)' 

Here,  for  the  tangential  component  of  the  contact  force,  we  are  content  with  the 
relationship  between  the  initial  increment  of  its  magnitude  T  and  that  of  the  tan¬ 
gential  displacement  s: 

T  =  Ks.  (7) 
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The  modulus  K  has  been  determined  by  Mindlin  [8]  for  the  initial  shearing  at  a 
contact  subsequent  to  the  application  of  the  normal  force  : 

K  =  ^^(3(1  -  (8) 

The  assumptions  are  that  there  is  Coulomb  friction  at  contacts  and  that  either  the 
relative  tangential  displacement  of  the  centers  of  contacting  spheres  is  very  small  or 
the  coefficient  of  friction  is  very  large.  In  either  case,  there  is  no  slip  at  contacts  and 
the  tangential  deformation  is  recoverable.  Following  Cauchy,  as  outlined  by  Love 
[7]  in  his  Note  B,  an  expression  for  the  average  stress  tensor  t  associated  with  a 
homogeneous  deformation  of  the  aggregate  may  be  obtained  by  considering  the  force 
transmitted  over  an  arbitrary  element  of  area  by  pairs  of  contacting  spheres  whose 
line  of  centers  is  cut  by  the  element.  For  random  assemblies,  the  stress  is  naturally 
expressed  in  terms  of  the  orientational  distribution  function  and  the  contact  force  : 

tij  =  )oijda  ,  (9) 

where  the  integration  is  over  all  solid  angle.  When  using  the  theory  to  interpret  the 
results  of  numerical  simulations  and  physical  experiments,  this  and  other  orienta¬ 
tional  averages  are  identified  with  volume  averages  over  sufficiently  large  regions  of 
nearly  homogeneous  strain  and  uniform  rotations. 

We  relate  the  displacement  u  of  a  contact  relative  to  the  center  of  a  sphere 
to  the  average  strain  of  the  aggregate  e  and  to  the  difference  between  the  average 
rotation  of  the  aggregate  w  and  the  average  rotation  of  the  spheres  a> : 

Ui  =  "b  (10) 

As  usual,  the  strain  tensor  is  symmetric  and  the  rotation  tensors  are  antisymmetric. 
As  we  shall  see,  it  is  important  to  make  the  distinction  between  the  two  average 
rotations  when  considering  the  static  theory  for  anisotropic  aggregates.  In  general, 
they  will  not  be  equal.  Schwartz,  Johnson,  and  Feng  [9]  and  Walton  [10]  made 
this  distinction  when  treating  the  dynamics  of  isotropic  assemblies,  but  found  that 
in  the  limit  of  long  wavelengths,  the  rotations  are  identical.  In  formulating  static 
theory  for  isotropic  aggregates,  Digby  [11]  and  Walton  [12]  tacitly  and  correctly 
assumed  that  the  rotations  were  the  same. 

Of  course,  even  with  the  incorporation  of  the  additional  degrees  of  freedom 
associated  with  the  average  rotations  of  the  particles,  the  kinematic  assumption 
(10)  may  be  far  from  correct.  For  example,  when  used  with  the  contact  force 
laws  (5)  and  (7),  it  requires  that  contacts  with  the  same  orientation  transmit  the 
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same  force.  Also,  when  the  average  strain  is  isotropic  and  compressive,  it  requires 
that  the  normal  components  of  all  contact  forces  must  be  equal  and  all  tangential 
components  must  vanish.  Numerical  simulations  described  by  Cundall,  Jenkins,  & 
Ishibashi  [13]  indicate  that  in  isotropic  aggregates  there  are  rather  large  variations 
in  the  forces  at  contacts  with  the  same  orientation,  even  when  the  mean  strain 
is  isotropic.  Also,  the  speed  of  the  shear  wave  through  a  compressed  isotropic 
aggregate  predicted  by  Digby  [11]  and  Walton  [12]  is  about  three  times  greater 
than  that  measured  in  numerical  simulations  and  experiments  [14].  However  in 
order  to  make  a  quantitative  evaluation  of  similar  shortcomings  of  the  kinematic 
assumption  (10)  in  anisotropic  aggregates,  we  must  first  obtain  the  predictions  of 
the  theory  that  is  based  upon  it. 

We  consider  a  transversely  isotropic  aggregate  that  has,  through  the  application 
of  the  appropriate  combination  of  isotropic  and  deviatoric  stress,  been  isotropically 
compressed  to  an  initial  average  volume  strain  of  magnitude  A.  In  this  initial  com¬ 
pression,  both  measures  of  the  average  rotation  are  assumed  to  be  zero.  Then,  we 
imagine  that  small  increments  of  homogeneous  strain  and  uniform  rotation  are  su¬ 
perposed  on  this.  In  this  event,  when  the  kinematic  assumption  is  used  to  calculate 
the  normal  and  tangential  components  of  the  contact  force,  we  obtain 


and 


P  =  MA^/^(A6ij  -  -eij)aiaj, 

Ti  =. 


(11) 


(12) 


These,  then,  are  used  with  the  orientational  distribution  function  (3)  in  the  expres¬ 
sion  (9)  for  the  stress  and  the  integrals  over  solid  angle  carried  out. 

Integrals  that  facilitate  this  calculation  axe 

hjki  =  JJ  aiOjOfea/da  =  ~{6ij6ki  +  6ik6ji  +  SuSjk)  (13) 

a 

and  the  more  elaborate 

// 

a 

The  result  of  the  integration  for  the  stress  is 


OliOljCXkOllOtpCiqdOC  —  ji^^iqljklp  "f"  ^jq^klpi  d"  ^kq^lpij  d"  ^Iqlpijk  d"  ^pqlijkl)-  (14) 


tij  =  d-  6ehihj]A 

127r  (1  —  I/) 


d- 

d- 


5  (2-1/) 

6ir  1/ 


1(5  2£)(€ij  ~h  tVij  ^ij')  d"  ^£hj{Cik  d"  Wik  i^ik^hk^ 


(15) 


[(7  4£)(efcfc(5tj  d"  -h  QEhkhlCkl^ij  d'  3£€kkhihj 


35  (2  -  u) 
d"  ^'^si.hiCjkhk  d"  hjCikhk]}  • 
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The  terms  in  the  first  line  of  equation  (15)  are  those  associated  with  the  initial 
isotropic  volume  strain.  Those  in  the  second  line  result  from  the  additional  incre¬ 
ments  of  strain  and  rotation  and  the  resistance  to  these  provided  by  the  component 
of  the  contact  force  that  is  tangent  to  the  surface  of  the  sphere.  The  terms  in  the 
third  and  fourth  lines  are  incremental  contributions  associated  with  both  compo¬ 
nents  of  the  contact  force.  If  the  contact  force  had  only  a  normal  component  with 
magnitude  given  by  (11),  then  the  stress  would  be  given  by  (15)  with  the  first  line 
unchanged,  the  second  deleted,  and  the  factor  i//(2  —  u)  multiplying  the  third  and 
fourth  fines  replaced  by  unity. 

When  body  couples  are  absent  and  the  deformation  is  homogeneous,  the  stress 
must  be  synunetric.  If  the  average  rotations  had  been  assumed  to  be  the  same  from 
the  outset,  symmetry  of  the  stress  could  be  achieved  by  discarding  the  antisymmet¬ 
ric  part  of  the  last  term  in  the  second  fine  of  (15).  The  resulting  expression  would 
be  symmetric,  but  incomplete.  The  correct  method  of  obtaining  a  symmetric  stress 
is  to  maice  the  distinction  between  the  rotations  and  to  determine  the  average  ro¬ 
tations  of  the  particles  so  that  the  symmetry  of  the  stress  is  insured.  For  example, 
symmetry  of  the  stress  given  by  (15)  requires  that 


2(5  -  2£){wij  -  ujij)  +  6e[(wijfc  -  uJik)hj  -  {wjk  -  u}jk)hi]hk 

-I-  Qe{eikhj  -  ejkhi)hk  =  0. 


(16) 


Solving  this  for  uf  yields 


+ 


3£ 

(5  -f  e) 


(Cikhj  Cjk^i)hk- 


(17) 


Thus,  given  the  strength  and  axis  of  the  anisotropy  and  the  average  strain  and 
rotation  of  the  aggregate,  equation  (17)  determines  the  average  rotation  of  the 
particles. 

Upon  employing  (17)  in  (15),  we  obtain  an  explicit  form  for  the  symmetric 
stress  given  in  terms  of  the  average  strain  of  the  aggregate  alone: 


6ehihj]A 

+  ^§5^1(5  -  2e)e„  +  +  h,ejt)hj, 

16^  /-  \  hkhlGklhihj 

(5  +  e) 


+ 

+ 


Gtt 


1/ 


[(7  ^E^Ckk^xj  d"  2eij)  "t"  GshkhlCkl^ij  d"  ScCfcfc /l.( /ij 


35  (2  -  u) 

\2£{hiejkhk  d-  h^Cikhk])  ■ 


(18) 
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This  is  the  form  of  the  stress  for  a  homogeneous  deformation  of  a  transversely 
isotropic  but  otherwise  random  array  that  has  been  first  isotropicedly  compressed 
and  then  subjected  to  a  small  strain.  Of  course  it  is  correct  only  in  so  far  as  the 
characterization  (3)  of  the  anisotropy,  the  contact  force  laws  (5)  and  (7),  and  the 
kinematic  assumption  (10)  are  correct. 


2  Inhomogeneous  Deformations 

When  the  deformation  of  the  aggregate  is  strongly  inhomogeneous,  the  calcu¬ 
lation  of  the  stress  is  slightly  more  complicated  and  a  corresponding  determination 
of  the  couple  stress  is  possible.  We  consider  an  infinitesimal  plane  element  of  area 
with  unit  normal  n  and  a  t5rpical  pair  of  contacting  particles  with  their  centers  on 
different  sides  of  the  plane.  The  vector  s  ~  (tol  is  directed  firom  the  center  of  the 
particle  in  the  region  pierced  by  — n  to  the  center  in  the  region  pierced  by  n.  The 
line  joining  the  centers  intersects  the  plane  at  the  point  r;  and  the  contact  is  located 
at  Tc  =  r  —  ^s,  where  — 1/2  <  ^  <  1/2.  The  location  of  the  area  element  is  assumed 
to  be  fixed  by  r  and  variations  in  the  position  of  the  point  of  intersection  over  the 
element  are  ignored.  The  number  density  n,  the  orientational  distribution  A,  and 
the  contact  force  F  are  all  taken  to  be  functions  of  Tc  and  a  .  Then 


A{rc,  -a )  =  A{r,  a  )  and  f^rc, -a  )  =  -F{rc,  a  )  (19) 

Agmn  the  determination  of  the  stress  parallels  that  of  Cauchy  as  outlined  by 
Love  [7].  However,  because  of  the  strong  inhomogeneity,  the  calculation  of  the  force 
transmitted  across  the  plane  by  all  pairs  of  particles  with  orientation  a  whose  line 
of  centers  is  cut  by  the  plane  involves  an  integration  over  The  result  is  most 
conveniently  expressed  in  terms  of 


fi{r  -  4s,a )  =  n{rc)A{rc,a)Fi{rc,a).  (20) 

Then,  at  r,  the  components  of  the  traction  r  ^u:e  given  in  terms  of  the  components 
of  n  by 

Tj  =  tijTlj,  (21) 

where 

Uj  f  //  /  “  )ajd^da  (22) 

with  the  double  integral  again  taken  over  all  solid  angle. 
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The  corresponding  expression  for  the  couple  stress  c  is  obtained  by  considering 
the  moment  of  the  contact  force  about  r: 


CiM  = -^Eijk  1 1  I  ^Sjfkir- ^8,OL)aid(da 
a 


(23) 


Then  the  components  of  the  couple  per  unit  area  x  at  r  are  given  in  terms  of  those 
of  c  and  n  by 

Xt  =  Ci/ni.  (24) 

The  stress  and  the  couple  stress  are  related  by  the  balance  of  moment.  In  the 
absence  of  body  couples,  this  is 


dcii 

dri 


^ijk^kj  —  0- 


(25) 


If  the  contact  forces  are  assumed  to  be  related  to  average  strains  and  rotations, 
then  simpler,  approximate  forms  for  the  stress  and  the  couple  stress  may  be  ob¬ 
tained.  These  axe  based  on  the  Taylor  series  expansion  of  the  contact  force  about 
r: 

/(r  -  fs)  =  /(r)  -  U  ■  V/(r)  +  .  V)  V(r).  (26) 

When  the  series  is  used  in  (22)  and  (23)  and  the  integrations  over  ^  are  carried 
out,  we  obtain 


Uj{r)  =  JJ  fi{r,a)ajda  ~  JJ  fi{r,a)ajakaidoi  (27) 

a  a 

and  ^  ^ 

Ciz(r)  = -— JJ  fk{r,a)ajaia^da,  (28) 

Q 

up  to  an  error  of  order  (cr,  L)^,  where  L  is  the  length  scale  over  which  the  average 
fields  vary. 

It  should  be  clear  from  (27)  and  (28)  that  the  couple  stress  and  the  antisym¬ 
metric  part  of  the  stress  depend  only  upon  the  tangential  component  of  the  contact 
force.  Also,  because  the  expansions  for  the  stress  and  couple  stress  have  been  car¬ 
ried  to  the  same  order,  there  is  an  additional  contribution  to  the  stress  that  involves 
the  second  spatial  derivatives.  It  may  be  that  this  term  is  small  in  comparison  to 
the  classical  part  of  the  stress,  however  the  contribution  of  its  antisymmetric  part 
to  the  balance  of  moment  (25)  is  exactly  the  same  order  as  the  divergence  of  the 
couple  stress  given  by  (28).  In  any  case,  it  is  the  diameter  of  the  particle  that 


provides  the  length  scale  in  terms  of  which  the  importance  of  the  new  terms  may 
be  evaluated. 

To  illustrate  the  calculation  of  the  couple  stress  and  new  parts  of  the  stress,  we 
suppose  that  the  aggregate  is  isotropic  and  the  behavior  of  a  contact  is  governed 
by  (4)  through  (8).  Then,  with  the  help  of  (13)  and  (14), 

^A{6ij6ki  +  SikSji  +  SiiSjk) 

a 


II 


fi{r,a)ajakQida  =  — 


(7  -  6u) 

(2-u) 


(^GijSkl  "I"  -|-  Cil6jk) 


V 

(2-1^) 


{,Gjk^il  “I"  GjlSik  Gkl^ij) 


-7 i^ijSki  +  +  ^iiSjk),  (29) 

where  O  =  w  —  CJ .  Using  this,  for  example,  in  (28)  an  exphcit  form  of  the  couple 
stress  may  be  obtained  : 


Cpk  -  [epjU^ijhi 

+Gpii(eik  +  fitfc)  +  £pki(Gii  +  fill)]-  (30) 

We  note  that,  as  a  consequence  of  the  assumptions  made  concerning  the  contact 
forces  and  the  kinematics,  the  gradient  of  the  difference  between  the  two  rotations 
appears  in  the  couple  stress,  often  accompanied  by  spatial  derivatives  of  the  strain. 
This  is  in  contrast  to  the  classical  Cosserat  theory  for  isotropic  materials  as  out- 
hned,  for  example,  by  Toupin  [15],  in  which  the  gradients  of  the  average  particle 
rotations  alone  appear.  Also,  in  this  classical  couple  stress  theory,  the  antisymmet¬ 
ric  part  of  the  stress  is  associated  with  the  difference  in  the  average  rotations  alone, 
not  with  the  spatial  gradients  of  the  difference  or  those  of  the  strain.  However,  the 
structure  of  the  theory  that  results  from  the  approach  outlined  above  does  depend 
rather  strongly  on  the  nature  of  the  contact.  For  example,  if  slip  rather  than  defor¬ 
mation  were  to  predominate  at  contacts,  then  particle  rotations  might  be  expected 
to  be  large  relative  to  the  classical  strains  and  rotations.  In  this  case,  the  resulting 
theory  for  the  inelastic  deformations  of  the  aggregate  is  likely  to  bear  a  stronger 
resemblance  to  the  Cosserat  theory.  In  any  event,  given  the  relation  of  the  contact 
forces  to  the  average  fields  of  interest  and  a  characterization  of  the  orientational 
distribution  of  contacts,  the  expressions  (27)  and  (28)  may  be  used  to  calculate  the 
stress  and  couple  stress.  Of  course,  the  relationship  between  the  results  of  this  cal¬ 
culation  and  what  is  measured  in  numerical  simulations  and  physical  experiments 
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will  depend  upon  how  closely  the  assumptions  upon  which  the  theory  is  based  are 
realized. 
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Chapter  Six 


Evolution  of  Elastic  Moduli 
in  a  Deforming  Granular  Assembly 

Peter  A.Cundall 


Abstract 


Numerical  experiments  are  performed  by  the  distinct  element 
method  on  a  sample  of  432  spheres  in  three  dimensions.  During  ax- 
isymmetric  loading  and  unloading  of  the  sample,  measurements  are 
made  of  incremental  elastic  moduli.  The  initial  elastic  shear  modulus 
compares  well  with  shear-wave  results  from  similar  physical  tests,  but 
the  reduction  in  modulus  with  shear  straining  is  greater  than  that  ob¬ 
served  physically.  The  evolution  of  confined  moduli  (corresponding  to 
pvwaves)  is  a  sensitive  measure  of  strain-induced  anisotropy. 


1  Background 

Chen  et  al.  (1988)  present  the  results  of  physical  experiments  on  samples  consisting 
of  a  mixture  of  two  sizes  of  glass  spheres.  The  experiments,  which  employed  a 
torsional  simple  shear  device,  produced  a  number  of  interesting  results  that  may 
be  interpreted  by  means  of  similar  numerical  experiments.  One  set  of  numerical 
experiments  was  performed  by  Cundall  (1988)  with  the  objective  of  reproducing 
the  stress /strain  behavior  in  shear  and  the  associated  volumetric  changes.  These 
experiments,  and  a  related  theoretical  study  by  Jenkins  (1988),  suggest  that  the 
physical  samples  contained  considerable  depositional  anisotropy. 

Similar  numerical  experiments  axe  reported  here,  but  the  emphasis  now  is  on 
the  variation  of  el£istic  wave  speeds  as  the  sample  is  loaded  and  unloaded. 

The  results  are  compared  to  the  physical  measurements,  and  are  related  to 
changes  in  sample  fabric. 


2  Numerical  Procedure 


The  distinct  element  method  (Cundall  &  Strack  1979a)  was  used  to  perform  the 
tests,  using  program  TRUBAL  (Cundall  &  Strack  1979b,  Strack  &  Cundall  1984)  to 
load  a  sample  contained  within  a  periodic  volume  (see  Cundall  1988).  A  432-sphere 
sample  was  prepared  with  the  following  characteristics  for  the  intact  material  (glass 
spheres)  and  for  the  assembly.  Note  that  the  elastic  and  frictional  properties  of 
glass  were  confirmed  by  direct  loading  on  individual  glass  beads  (Ishibashi  1989). 


Shear  modulus  of  glass 
Poisson’s  ratio  of  glass 
FYiction  coefficient 
Number  of  particles: 

Isotropic  stress 
Initial  porosity 
Initial  coordination  number 


G  =  2.9  X  10^°  Pa 
u  =  0.2 
fi  =  0.3 

40  of  0.1825mm  rad. 
392  of  0.1075mm  rad, 
<To  =  1.38  X  10^  Pa 
no  =  0.368 
Ko  =  5.37 


Hertz  theory  applies  in  the  normal  contact  direction.  In  the  shear  direction,  the 
stiffness  depends  on  normal  force  but  not  on  shear  displacement  (see  Cundall  1988). 
In  order  to  achieve,  simultaneously,  a  given  isotropic  stress  and  a  given  porosity, 
some  adjustment  of  the  friction  coefficient  is  necessary  during  compaction.  The 
following  sequence  of  steps  was  used  to  create  the  initial  isotropic  sample. 

1.  Random  placement  of  all  particles  within  a  cubic  space  of  19  mm^;  no  particle 
touches  any  other. 

2.  Isotropic  shrinking  of  the  periodic  volume  to  a  porosity  of  n  =  0.367,  under  a 
friction  coefficient  of  =  0.01. 

3.  Isotropic  compaction  or  expansion  of  the  periodic  volume  under  n  =  0.05,  until 
a  mean  stress  of  Oo  =  1-38  x  10^  Pa  is  achieved. 

4.  Setting  of  friction  coefficient  /i  to  0.15,  and  further  isotropic  volume  adjustments 
eis  necessary  to  keep  the  meein  stress  at  1.38  x  10®  Pa.  The  final  porosity  stabilized 
at  n  =  0.368.  The  friction  coefficient  was  then  set  to  0.3. 

For  comparison,  the  physical  tests  (Chen  et  al.  1988)  were  done  under  conditions 
of  (To  =  138  KPa,  n  =  0.367. 

Compression/extension  shear  tests  are  performed  by  applying  strain  rates  to 
the  periodic  space  as  follows 
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where  'ym  is  the  maximum  shear  strain  plotted 
in  the  figures  that  follow.  Constant  isotropic 
stress  is  maintained  with  a  numerical  servo- 
control  that  adjusts  e^. 

3  Results 

Figure  1  provides  three  orthogonal  views 
of  the  sample  after  volumetric  compaction:  the 
intensity  of  shading  is  related  to  the  mean 
stress  carried  by  each  particle.  White  (tm- 
shaded)  particles  carry  no  load  -  i.e.  they  are 
“floating.”  It  is  clear  that  there  is  a  wide  varia¬ 
tion  in  the  average  load  carried  by  the  particles, 
even  though  the  compaction  was  carried  out  at 
a  very  low  firiction  coefficient. 

A  histogram  of  normal  forces  is  presented 
in  Figure  2  for  the  same  starting  state. 

Figure  3  records  the  stress/strain  response 
of  the  sample  for  axisymmetric  loading  and  im- 
loading  according  to  equations  (1). 

Figure  4  shows  the  corresponding  volu¬ 
metric  response,  and  Figiire  5  the  evolution  of 
coordination  number,  K,  with  strain. 

The  maximum  applied  strain  was  about  0.45%,  to  correspond  with  the  physical 
tests  performed  by  Chen  et  al.  (1988)  and  Ishibashi  et  al.  (1988). 

In  the  physical  tests  of  Chen  et  al.  (1988),  long  wavelength  shear  waves  of 
low  amplitude  were  propagated  through  the  samples  in  order  to  assess  the  chang¬ 
ing  elastic  shear  modulus.  For  low  amplitude  oscillations,  contact  sliding  does  not 
occur — sliding  stops  at  the  instant  of  unloading;  upon  reloading,  sliding  only  re¬ 
sumes  again  at  the  original  point  of  loading. 


Figure  1:  Three  orthogonal  views 
of  the  numerical  assembly  at  its 
initial  state;  the  intensity  of  shad¬ 
ing  is  proportional  to  the  mean 
stress  carried  by  each  sphere. 
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Normal  force  (in  30  ranges) 


Figure  2:  Histogram  of  contact  normal 
forces  at  the  initial  state  (A). 


forces  at  the  initial  state  (A). 


Figure  4:  Volume  changes  corresponding  to  Figure  3. 

Such  “probes”  are  simulated  numerically  by  making  small  strain  increments, 
Aei2,  at  infinite  friction  to  prevent  sUding.  The  probes  were  done  at  several  points 
along  the  loading/unloading  path,  and  corresponding  shear  modulus  values  plotted 
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Maximum  shear  strain:  % 

Figure  5:  Evolution  of  coordination  number  during  loading  and  unloading  test  of 
Figure  3. 


0.0  0.1  0.2  0.3  0.4  0.5 

Maximum  shear  strain:  % 

Figure  6:  Evolution  of  elastic  shear  modulus  (in  ei2  direction)  during  the  test  of 
Figure  3. 


The  imtiai  modulus  (at  zero  shear  strain)  is  127  MPa,  which  compares  with  161 
MPa  derived  from  propagating  shear  waves  in  the  physical  tests  (Chen  et  al.).  The 
reduction  in  numerical  shear  modulus  at  0.45%  strain  is  32%,  while  the  physical 
tests  give  a  reduction  of  around  20%  for  the  same  shear  strain. 

It  is  known  that  a  granular  material  becomes  elastically  anisotropic  when 


loaded  in  shear.  Partly,  this  is  due  to  preferential  contact  loss  in  the  minor  prin¬ 
cipal  strain  direction:  i.e.  a  geometric  fabric  develops.  But  when  contacts  obey 
the  Hertzian  law,  they  become  stiffer  in  the  major  direction  and  softer  in  the  mi¬ 
nor  direction:  this  dependence  on  the  angular  force  distribution  may  be  termed 
a  kinetic  fabric.  Both  geometric  and  kinetic  effects  occur  in  the  numerical  tests, 
but  there  does  not  appear  to  be  an  accepted  way  to  measure  the  combined  fabric. 
We  can  quantify  the  anisotropy  by  making  small-strain  probes  that  correspond  to 
propagating  compressional  waves  in  three  orthogonal  directions;  recall  that  en  is 
in  the  major  principal  direction  of  the  axisymmetric  sample.  Table  1  records  the 
modulus  values  at  various  states  of  the  loading/umoading  test — refer  to  Figure  3. 
The  probes  are  applied  as  strain  increments  in  the  given  directions;  all  other  strain 
components  are  zero. 


State 

(11) 

(22) 

(33) 

A 

341.1 

357.6 

356.9 

B 

366.0 

204.5 

173.1 

C 

310.4 

279.7 

248.7 

D 

157.5 

289.8 

307.1 

Table  1:  Confined  modulus  values  in  MPa  for  probes  in  three  orthogonal  directions. 


Figures  7  and  8  record  the  distributions  of  contact  normal  forces  at  states  B 
and  C,  respectively. 
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(mean)  (max) 

Normal  force  (in  30  ranges) 


Figure  7:  Histogram  of  contact  normal  forces  at  state  B  (maximum  shear  strain). 


Figure  8:  Histogram  of  contact  normal  forces  at  state  C  (zero  shear  stress,  after 
imloading). 

4  Conclusions 

The  agreement  in  initial  elastic  shear  moduli  between  physical  and  numerical  tests 
is  considered  to  be  quite  good,  bearing  in  mind  that  the  numerical  simulations  were 
bgtsed  on  properties  of  intact  glass,  and  incorporated  no  arbitrary  factors.  However, 
the  numerical  tests  show  a  50%  greater  reduction  in  modulus  with  straining  than 
the  physical  tests.  Table  1  suggests  that  p-wave  measurements  in  tests  on  granular 
material  may  be  used  as  a  sensitive  measure  of  anisotropy,  and  how  it  evolves  with 
shear  strain.  The  ratio  between  orthogonal  modulus  values  is  almost  2:1  at  0.45% 
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shear  strain. 

The  histograms  of  normal  force  show  that  force  distributions  are  heavily  biased 
towards  the  low  values;  theories  that  assume  symmetrical  distributions  about  the 
mean  may  be  seriously  in  error. 

The  simulations  were  all  done  on  an  80386-based  micro-computer,  8md  each 
stage  (e.g.  compaction;  loading)  takes  two  to  four  hours  to  compute;  each  probe 
takes  25  minutes. 
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Chapter  Seven 


Measurement  of  Local  Strainrate 
in  a  TRUBAL  Sample 


Peter  A.  Cundall 


Abstract 


The  numerical  simulation  TRUBAL  is  employed  together  with 
a  least-squares  procedure  to  determine  the  local  strain  tensor  expe¬ 
rienced  by  a  particle.  The  average  of  the  straing  tensor  determined 
from  the  local  strain  tensors  is  compared  to  the  overall  strain  tensor 
the  difference  between  the  two  is  calculated. 


1  Discussion 

Most  theoretical  models  of  a  granular  material  contain  the  assumption  that  a 
representative  particle  is  embedded  in  a  uniform  strain  field,  in  order  to  compute 
the  relative  displacements  at  particle  contacts.  It  is  instructive  to  use  TRUBAL  to 
measure  the  local  strain  tensor  experienced  by  each  particle,  and  develop  statistics 
that  express  the  spacial  fluctuations  in  strainrate. 

The  measurement  of  local  strain  is  not  as  simple  as  the  measurement  of  local 
stress.  To  determine  the  average  stress  tensor  of  a  particle  we  use  an  expression  that 
contains  the  discrete  forces  acting  at  the  particle’s  contacts;  the  forces  elsewhere 
are  zero.  It  is  not  correct  to  use  the  displacements  in  a  similar  way  to  find  strains, 
because  the  displacements  in  the  voids  are  not  zero.  The  method  adopted  here 
is  to  find  the  strainrate  tensor  that  minimizes  the  error  between  the  predicted 
displacements  and  the  measured  displacements.  A  least-squares  procedure  is  used 


to  derive  the  tensor,  given  a  set  of  displacement  vectors  and  corresponding  position 
vectors.  The  general  procedure  is  described  first,  and  tested  on  manufactured  data, 
consisting  of  displacements  at  random  points  derived  from  random  components  of  a 
displacement  gradient  tensor.  It  should  be  noted  that,  in  general,  the  displacement 
vector  components  do  not  sum  to  zero  —  i.e.  there  is  a  rigid-body  movement.  The 
least-squares  procedure  also  determines  this. 


1  General  Least- Squares  Procedure 

We  wish  to  find  the  displacement  gradient  tensor,  Cij,  and  the  rigid-body  dis¬ 
placement,  u°,  that  represent  the  best  fit  to  a  set  of  n  measured  displacement  values, 
up,  where  m  =  1,  n.  The  displacement  vector  Up  is  located  at  a  coordinate  of  x^. 
The  predicted  displacements,  up  at  the  n  points  are: 


=  P 


I  A  *  ^ 


(3) 


A  measure  of  the  error  is  the  sum  of  the  squares  of  the  deviations  between  predicted 
and  measured  displacements: 


m=l  i=l 

The  condition  for  minimum  z  is  that 


— —  =  0  and 

dcij 


(4) 

(5) 


Substituting  (3)  into  (4)  and  differentiating,  a  set  of  twelve  equations  is  obtained 
as  follows,  where  i  =  1,3: 


e,,Y.^riT  +  e.2  5]  +  e.3  +  <  E  =  E  (®) 

m  m  m  mm 

e,i  Y,  +  ei2  Y  +  ea  Y  +  <  E  =  E  (’’) 

m  TTi  m  m  m 

Cii  ^  xpxp  +  ei2  xpxp  +  ei3  x^xp  +  <  =  51  Upxp  (8) 

TTI  m  m  m  m 

ei,53ir  +  ei2E^"  +  «-3E^"  +  <  =  E^r  (9) 

m  m  m  m 
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Since  the  coefficients  of  Cij  and  u°  on  the  left-hand  sides  are  independent  of  we 
need  only  evaluate  sixteen  coefficients,  and  then  apply  a  solution  procedure  three 
times,  with  different  right-hand  sides.  The  LU-decomposition  routine  LUDCMP  given 
in  Press,  et  al.  (1986)  is  applied  to  the  4x4  coefficient  matrix,  and  the  accompanying 
back-substitution  routine  LUBKSB  used  three  times;  in  this  way,  all  twelve  unknowns 
(nine  strains  and  three  displacements)  are  obtained. 

The  scheme  is  tested  by  starting  with  random  strain  and  rigid-body  compo¬ 
nents.  These  are  used  to  generate  displacement  vectors  at  n  random  locations, 
using  (3).  The  solution  process  described  above  then  produces  computed  strain 
and  rigid-body  components:  these  are  compared  to  the  assumed  starting  values. 
It  is  found  that  the  computed  and  assumed  values  are  very  close  for  n  >  4.  The 
program  TEST,  listed  in  Appendix,  performs  this  comparison. 

2  Least-Squares  Procedure  Applied  to  TRUBAL  Data 

m  The  above  procedure  is  applied  to  sets  of  velocity  vectors  derived  from  a 
sphere  assembly  undergoing  continuous  straining.  For  each  sphere  (the  “target” 
sphere)  in  the  assembly,  all  contacting  spheres  are  identified.  Two  separate  ways  of 
providing  velocity  data  are  used  (yielding  two  distinct  measures  of  local  strainrate): 

1.  The  translational  velocities  of  contacting  spheres  relative  to  the  target 
sphere  are  taken;  there  is  no  contribution  from  the  spins  of  any  particle. 
The  coordinates  inserted  into  equations  6-9  are  the  relative  centroid 
coordinates. 

2.  At  each  contact  on  the  target  sphere,  the  relative  velocity  vector  between 
two  coincident  points  on  the  two  spheres  is  used  in  the  equations.  This 
strain  measure  is  affected  by  the  spins  of  the  target  particle  and  all 
contacting  particles.  The  coordinates  of  the  contact  points  are  used, 
rather  than  the  centroids. 

The  procedure  outlined  above  is  implemented  into  subroutine  JSTAT,  and  called 
from  TRUBAL.  It  was  checked  with  a  number  of  two-  and  five-particle  simulations,  in 
which  the  spheres  are  given  fixed  v  elocities.  Various  combinations  of  rotations,  spins 
and  translations  produced  strainrates  and  rigid-body  motions  that  were  compared 
with  predicted  values.  At  this  level,  the  procedure  appeared  to  work  as  intended. 

3  Preliminary  Rresults  Prom  TRUBAL  Simulations 

Some  preliminary  tests  have  been  done  with  the  54-sphere  assembly  used  in 
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previous  tests  (Cundall,  1988).  These  tests  use  Hertzian  contact  laws.  An  unex¬ 
pected  result  is  that  significant  relative  rigid-body  motions  are  observed  —  i.f;.  the 
target  particle  appears  to  have  a  continuous  translational  velocity  relative  to  its 
contacting  neighbors.  At  first  sight  this  does  not  seem  to  be  reasonable.  However 
the  combination  of  anisotropic  contact  distribution  and  a  nonlinear  contact  law 
means  that  overall  tangent  stiffnesses  vary  with  angle  around  a  sphere.  If,  say,  one 
principal  stress  is  increased,  the  relative  contact  displacements  on  opposite  sides  of 
a  sphere  must  be  different,  if  equilibrium  is  to  be  maintained.  Hence  the  target 
sphere  moves  relative  to  its  neighbors. 

The  following  results  are  for  four  stages  in  an  axisymmetric  compression  test, 
starting  with  the  isotropic  state  produced  by  data  file  SRI. DAT.  The  sample  is 
compressed  in  the  cn  direction,  using  a  servo-control  to  keep  the  mean  stress  cto 
constant  at  138  KPa.  The  shear  stress  ratio,  t/<To,  is  shown  for  each  stage,  where 

(T22  +  0~33  1 

2  J 


The  applied  strainrate  components  (of  the  periodic  space)  are  shown  first,  in  the 
form 


eii 

612 

613 

621 

622 

623 

631 

632 

633 

The  strainrates  represent  averages  of  computed  values  for  all  particles  that  have 
four  or  more  contacts:  the  actual  number  of  particles  is  shown.  The  two  measured 
strainrate  tensors  are  then  given,  in  the  same  form,  followed  by  the  relative  rigid- 
body  velocities  for  the  two  cases  discussed  previously.  The  velocities  are  given  in 
the  form:  (1*1,112,^3). 


State  at  cycle  29050;  tcto  =  0.31. _ 

Applied  strainrates  .  .  . 

-1.304e-07  0  0 

0  7.956e-08  0 

0  0  7.956e-08 

Strainrate  measure  1  ... 
-1.317E-07  -2.702E-09  -1.365E-08 
-l.lOOE-07 

7.851E-09  7.306E-08  -1.185E-09 

-2.181E-08 

-9.145E-09  -3.549E-09  7.914E-08 

5.862E-08 


(4o  particles  considered) 
Strainrate  measure  2  . . . 
-1.220E-07  -1.166E-08 

7.623E-09  6.717E-08 

7.862E-08  7.423E-09 
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Relative  velocities  . . . 


1.606E-06  1.273E-06  1.306E-06  1.265E-06  9.717E-07 


9.694E-07 


State  at  cycle  30050;  t/cto  =  0.42. 

Applied  strainrates  . . . 

-1.296e-07  0 

0 

0  8.042e-08 

0 

0  0 

8.042e-08 

(47  pzurticles  considered) 

Strainrate  measure 

1  . . . 

Strainrate  measure  2  ... 

-1.385E-07  4.289E-08 

-1.005E-07 

-7.000E-08 

-1.153E-07  5.952E-09 

1.570E-08  8.747E-08 

-4.322E-08 

1.317E-08 

1.233E-08  5.250E-08 

1.959E-08  -4.749E-08 

4.659E-08 

Relative  velocities  . . . 

1.065E-07 

6.215E-08  3.052E-08 

1.898E-06  2.127E-06 

1 . 248E-06 

1.998E-06 

1.620E-06  1.124E-06 

State  at  cycle  31050;  r/<To  =  0.47. 


Applied  strainrates  . . . 

-1.263e-07  0 

0 

• 

0  8.370e-08 

0  0 

0 

8.370e-08 

(45  particles  considered) 

Strainrate  measure 

1  . .  . 

Strainrate  measure  2  . . . 

-1.374E-07  1.033E-08 

-3.049E-08 

-1.072E-07  3.644E-09 

-3.847E-08 

-1.388E-08  6.290E-08 

1 . 228E-08 

-l,917E-08  3.997E-08 

1.332E-08 

3.278E-08  2.573E-08 

8.349E-08 

1.102E-07  6.303E-08 

6.426E-08 

Relative  velocities  . .  . 

2.246E-06  2.606E-06 

2.842E-06 

2 . 233E-06  1 . 502E-06 

1.707E-06 

State  at  cycle  32050;  t/(To  —  0.49. _ 

Applied  strainrates  . . . 

-1.259e-07  0  0 

0  8.408e-08  0 
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0 


0  8.408e-08 


Strainrate  measure  1  ... 
-1.349E-07  2.004E-08  -3.243E-08 

-1.774E-09 

-3.121E-08  6.634E-08  1.840E-08 

-1.067E-08 

2.262E-08  5.704E-08  8.106E-08 

1 . 072E-07 

Relative  velocities  . . . 

2.732E-06  3.337E-06  3.007E-06 

2.401E-06 


(44  particles  considered) 
Straijuxate  measure  2  .  .  . 
-2.061E-07  4.816E-08 

-2.584E-08  6.380E-08 

7.896E-08  9.321E-08 

2.984E-06  1.677E-06 


The  main  observation  concerning  the  results  is  that  the  ensemble  averages  do 
not  correspond  very  well  with  the  appUed  strainrates;  furthermore,  there  are  large 
asymmetric  components.  It  seems  pointless  to  try  to  study  fluctuations  in  strain 
(i.e.  deviations  from  the  average),  when  the  averages  are  themselves  exhibit  such 
noisy  behavior.  The  noise  may  come  from  an  error  in  the  procedure  or  coding,  or 
it  may  be  real.  In  either  case,  further  work  is  needed  to  gain  understanding. 
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APPENDIX:  listing  of  TEST  program  for  strain  calculation 

program  test 
c 

c  Test  of  least-squares  recovery  of  strainrate  tensors  from 

velocities 

c 

logical  flag 
parameter  (ns=4,np=10) 

dimension  a(ns,ns) ,b(ns) ,indx(ns) ,etarg(3,3) ,ecomp(3,3) 
dimension  x(3,np) ,v(3,np) ,rbv(3) 

10  write  (*,1001) 
read  (*,*)  num 
if  (num  .eq.  0)  stop 

c -  random  rigid-body  motions  - 

rbv(l)  =  (uranO  -  0.5)  ♦  le-5 

rbv(2)  =  (uranO  -  0.5)  *  le-5 

rbv(3)  =  (ureinO  -  0.5)  *  le-5 

write  (*,1006) 

c -  generate  target  strain  tensor  - 

do  20  i  =  1,3 

do  15  j  =  1,3 

etarg(i,j)  =  le-6  *  (uranO  -  0.5) 

15  continue 

write  (*,1002)  (etargd,  j)  ,  j=l  ,3)  ,rbv(i) 

20  continue 

write  (*,1004) 

c -  generate  measurement  points  - 

do  40  n  =  1 , num 
do  25  i  =  1,3 

x(i,n)  =  10.0  *  uranO 
25  continue 

do  35  i  =  1,3 
v(i,n)  =  rbv(i) 
do  30  j  =  1,3 

v(i,n)  =  v(i,n)  +  etargd, j)  *  x(j,n) 

30  continue 
35  continue 

write  (*,1003)  n, (x(i,n) ,i=l,3) ,  (v(i ,n) , i=l ,3) 
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40  continue 

flag  =  .false. 

c -  assemble  A  matrix  - 

do  60  j  =  1,3 

do  50  i  =  1,3 

a(i,j)  =  0.0 
do  45  n  =  l,num 

a(i,j)  =  a(i,j)  +  x(i,n)  ♦ 
45  continue 
50  continue 
60  continue 

do  68  k  =  1,3 
a(4,k)  =  0.0 
a(k,4)  =  0.0 
do  65  n  =  l,num 

a(4,k)  =  a(4,k)  +  x(k,n) 
a(k,4)  =  a(k,4)  +  x(k,n) 

65  continue 
68  continue 

a (4, 4)  =  float (num) 

c -  LU  decomposition  - 

call  ludcmp  (a,ns,ns,indx,d,flag) 
if  (flag)  then 
write  (*,1000) 
goto  500 
else 

write  (*,1007) 

c -  assemble  b  vector  for  each  axis  — 

err  =  0.0 
do  85  irow  =  1,3 
do  75  i  =  l,ns 
b(i)  =  0.0 
do  70  n  =  l,num 

if  (i  .eq.  4)  then 

b(i)  =  b(i)  +  v(irow,n) 
else 

b(i)  =  b(i)  +  v(irow,n) 
endif 


x(j  ,n) 


*  x(i,n) 


75 


70  continue 

75  continue 

call  lubksb  (a,ns,ns,indx,b) 
write  (♦,1005)  (b(i) ,i=l ,ns) 

85  continue 
endif 

500  goto  10 

1000  format  (’  *♦♦  singular  matrix’) 

1001  format  (’  Number  of  points  (0  to  quit)  ?’$) 

1002  format  (4x,lp,3ell .3,3x,ell .3) 

1003  format  (lx,i3,lp,3ell.3,5x,3ell.3) 

1004  format  (’  Random  locations;  computed  velocities  ...’) 

1005  format  (4x, lp,3ell .3,3x,ell .3) 

1006  format  (’  Target  strains  k  rigid-body  velocities  ...’) 

1007  format  (’  Generated  strains  k  rigid-body  velocities  ...’) 
end 

function  uran  () 
c 

c  Uniform  random  niunber  generator,  according  to  Wichmein  k  Hill, 
c  Byte,  March  1987. 
c  Seed  is  set  in  data  statement, 
c 

save  ix,iy,iz 

data  ix,iy,iz  /I , 10000,3000/ 
j  =  ix  /  177 
k  =  ix  -  177  ♦  j 
ix  =  171  ♦  k  -  2  ♦  j 
if  (ix  .It.  0)  ix  =  ix  +  30269 
j  =  iy  /  176 
k  =  iy  -  176  ♦  j 
iy  =  172  ♦  k  -  35  ♦  j 
if  (iy  .It.  0)  iy  =  iy  +  30307 
j  =  iz  /  178 
k  =  iz  -  178  ♦  j 
iz  =  170  *  k  -  63  *  j 
if  (iz  .It.  0)  iz  =  iz  +  30323 
temp  =  float (ix)/30269.0  +  float (iy)/30307.0  + 
float (iz)/30323.0 
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Chapter  Eight 


Mean-field  Inelastic  Behavior  of 
Random  Arrays  of  Identical  Spheres 


James  T.  Jenkins  Otto  D.  L.  Strack 


Abstract 


We  consider  a  random  array  of  identical  spheres  that  interact 
through  noncentral  contact  forces.  We  assume  that  the  displacement 
of  a  contact  relative  to  a  center  may  be  calculated  from  the  average 
strain  of  the  aggregate.  The  normal  component  of  the  contact  force  is 
assumed  to  be  Hertzian  and  the  tangential  component  is  assumed  to  be 
linearly  elastic  until  frictional  sliding  occurs.  We  consider  the  response 
of  the  material  in  triaxial  compression.  For  monotone  deformations,  we 
calculate  the  evolution  of  the  contact  distribution,  the  volume  change, 
the  stress-strain  response,  the  plastic  strain,  and  the  strain  hardening. 


1  Introduction 

We  are  interested  in  the  average  small  strain  response  of  random  granular  as¬ 
semblies  that  interact  through  contacts  that  may  slide.  This  interest  originated  in 
attempts  to  explain  the  stress-strain  behavior,  volume  changes,  and  wave  propaga¬ 
tion  observed  in  deformations  of  a  granular  material  carried  out  at  a  fixed  pressure. 
These  experiments  were  done  on  hollow  cyUndrical  samples  of  essentially  identical 
glass  spheres  confined  within  a  membrane  in  a  true  triaxial/torsion  device  (Chen, 
Ishibashi,  &  Jenkins,  198).  Consequently,  we  focus  our  attention  on  an  idealized 
material  consisting  of  identical  spherical  grains  and  assume  that  the  location  of 
their  centers  is  random. 

In  triaxial  compression  and  extension,  we  assume  that  the  deformation  of  a 


sphere  in  the  neighborhood  of  a  contact  can  be  expressed  in  terms  of  the  average 
strain  of  the  aggregate.  We  are  aware  that,  because  of  irregularities  in  the  ar¬ 
rangement  of  the  particles,  differences  in  the  nvunber  of  contacts  per  particle,  and 
variations  in  the  magnitude  of  the  contact  forces,  this  mean  field  assumption  is  only 
approximately  true.  However,  we  feel  that  it  is  worthwhile  to  derive  a  complete  a 
theory  based  on  it  in  order  to  determine  its  structure  and  to  obtain  its  predictions. 
The  structure  of  more  complicated  theories  can  be  expected  to  contain  the  elements 
of  the  simplest  theory,  and  the  predictions  of  the  simplest  theory  can  be  obtained 
analytically. 

2  Theory 

We  suppose  that  the  diameter  of  a  sphere  is  D  and  that  there  are  n  of  them 
per  unit  volume.  We  take  a  to  be  the  unit  vector  from  the  center  of  a  sphere  to 
a  contact  on  its  surface.  The  rectangular  Cartesian  components  of  the  unit  vector 
a  axe  (cos  0  sin  0,  sin  0  sin  0,  cos  6),  where  9  is  the  polar  angle  from  the  axis  of 
symmetry. 


3  Contact  Displacement 


As  mentioned  above,  we  assume  that  the  displacement  u  of  a  contact  point 
relative  to  the  center  of  its  sphere  is  given  in  terms  of  the  average  strain  e  of  the 
aggregat'^  and  the  vector  from  the  center  of  the  sphere  to  the  contact  by 

Ui  —  2  CijCkj.  (1) 

In  more  complicated  deformations  involving  average  rotations,  the  corresponding 
mean  field  assumption  involves  both  the  average  rotation  of  the  aggregate  and  the 
average  spin  of  the  spheres  (Jenkins,  1991). 

In  triaxial  compression  or  extension,  the  volume  strain  A  (taken  positive  for  a 
decrease  in  volume)  is 

A  =  —(633  +  2eii)  (2) 

and  the  shear  strain  7  (half  the  mzucimum  shear  strain)  is 


7  =  “2(^33  -  eii). 

The  normal  displacement  ^  of  a  contact  point  toweird  a  center  is,  from  (1), 

^  =  -^(A  —  27  -f  67cos^  6). 


(3) 

(4) 
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The  corresponding  tangential  displacement  s  is 


s  =  £)7sin0cos0ea,  (5) 

where  ee  is  the  unit  vector  in  the  direction  of  increasing  9. 

4  Contact  Force 

The  force  F(a )  exerted  by  the  sphere  at  a  contact  with  orientation  a  has 
components  parallel  and  perpendicular  to  a : 

Fi=Pai-Ti,  (6) 

where  T-  a  =0.  In  a  homogeneous  deformation,  F  is  independent  of  position. 
Then,  because  the  force  exerted  on  the  sphere  at  a  contact  is  equal  and  opposite  to 
the  force  exerted  by  the  sphere,  F(-a  )  =  -J^(a ).  The  magnitudes  of  the  normal 
and  tangential  components  depend  upon  the  displacement  of  a  contact  relative  to 
the  center  and  are,  in  general,  functions  of  a . 

Following  Hertz,  the  normal  component  P  is  related  to  the  normal  component 
of  the  displacement  (Love,  1927): 

P  =  (7) 

where  M  is  given  in  terms  of  the  shear  modulus  G  and  Poisson’s  ratio  of  the 
material  of  the  spheres  by 

2  GD^ 

M  =  - ■=— - r.  (8) 

9n/3  (1  - 1") 

Here,  for  the  tangential  component  of  the  contact  force,  we  employ  a  bilinear 
relationship  that  incorporates  elastic  displacement  and  frictional  sliding.  Its  form 
depends  on  the  magnitude  T  of  the  tangential  component  relative  to  the  product  of 
the  coefficient  of  friction  /x  between  the  spheres  and  the  normal  force  P.  Provided 
that  T  <  fiP,  T  is  related  to  the  magnitude  s  of  the  tangential  displacement  through 

T  =  Ks,  (9) 


where  the  modulus  K  is  that  determined  by  Mindlin  (1949)  for  the  initial  shearing 
at  a  contact  subsequent  to  the  application  of  the  normal  force: 

05/3/^2/3 

K  =  |3(1  -  y)DP]''K  (10) 
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Otherwise, 


T  =  nP. 


(11) 


Equations  (9)  through  (11)  are  an  integrated  bilinear  approximation  to  the  more 
elaborate  incremental  inelastic  behavior  discussed  by  Mindlin  and  Deresiewicz 
(1953). 


5  Contact  Deletion  and  Sliding 

The  results  concerning  contact  deletion  and  sliding  in  triaxial  compression  (7  > 
0)  are  most  compactly  expressed  in  terms  of 


6  =  A  —  27  =  — 3eii  and  c  =  67. 


(12) 


Also,  it  is  convenient  to  restrict  our  attention  to  the  range  0  <  ^  <  7r/2. 

At  the  onset  of  sliding,  the  tangential  displacement  is  completely  elastic. 
With  the  help  of  (9),  (11),  (7),  (8)  and  (10),  its  magnitude  s®  is  given  by 

=  T/K  =  ^lP/K  =  /i<5,  (13) 


(14) 


where 

^  '  3  1  -  !/■ 

When  (5)  and  (4)  are  used  in  (13),  the  condition  for  sliding  to  begin  may  be  ex¬ 
pressed  in  terms  of  6,  c,  and  9  as 


c  sin  0  cos  0  =  fi{b  +  c  cos^  9). 


(15) 


Sliding  does  not  occur  while  the  function  g{9),  defined  as 

g{9)  =  s  —  p,6  =  ^[c(sin  29  —  p,  cos  29  —  p)  —  2pb]  (16) 

is  negative.  The  first  zero  for  g{9)  occurs  in  the  interior  of  the  interval  at  the  value 
9c  for  which  dg/d9  =  0.  We  have,  from  (16),  that 

^  =  ^  (cos  29  c  +  P  sin  29c)  =0  (17) 

da  6 

or 

cot  20c  =  — A  (18) 

Note  that  9c  is  independent  of  b  and  c  and  greater  than  7r/4. 
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Substitution  of  9c  for  6  in  g{dc)  =  0  yields  the  relation  between  h  and  c  at  the 
onset  of  sliding: 

OZ) 

(19) 


COS  20c 

c  =  -26- 


1  -I-  cos  26  c 


With  (18),  (16)  may  be  written  as 


M  ^  cot  20c 


{•[ 


cos  2(0  —  0c) 
cos  20c 


-1-1 


-I-  26 


}■ 


In  further  shearing, 


c>  -26 


cos  20c 
1  -I-  cos  20c 


>  0. 


(20) 


(21) 


Sliding  then  occurs  over  a  range  of  angles  6m  <  9  <  6m,  where, 
and  6m  axe  obtained  as  the  roots  of  ^(0)  =  0: 

when  b  >  0,  6m 

cos  2(0  -  0c)  =  —  (l  -1-  2-)  cos  20c. 

(22) 

The  roots  are 

6m  =  6c  —  \  arccos  [(l  +  2-)  cos  20c] 

(23) 

and 

^  arccos  [(l  +  2-)  cos  20c] . 
c 

(24) 

Contact  is  effectively  lost  at  those  contacts  at  which  the  normal  component  of 
the  contact  force  has  relaxed  to  zero  or,  equivalently,  when  6,  given  by  (4),  vanishes. 
In  triaxial  compression,  contact  is  first  lost  at  0  =  7r/2  when  6  =  0.  In  order  to 
characterize  further  loss  of  contact,  we  define  an  angle  0i  by 

r  7r/2,  when  6  >  0 

6i  =  's  /— — 

[  arccos  y/ —b/c,  when  6  <  0. 

(25) 

Then,  when  6  <  0,  contact  is  lost  over  the  range  in  angle  6\  <6  <-k/2  and  sliding 
takes  place  over  the  range  6m  <0  <  6m,  with  6m  given  by  (24)  and  6m  =  0\. 

When  sliding  occurs,  the  magnitude  of  the  part  of  the  tangential  displace¬ 

ment  associated  with  it  is  equal  to  the  function  ^(0): 

=  g{6)  =  [c  sin  0  cos  0  —  p{b  +  c  cos^  0)] , 

6 

(26) 

where  6m  ^  0  <  6m- 

For  glass,  typical  values  of  /z  and  u  are  0.3  and  0.2,  respectively,  so  p,  =  0.225. 
This  value  of  p  is  used  when  displaying  the  results  of  the  calculations. 
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6  Contact  Distribution 


Let  duj  =  sin  6d6d4>  be  the  element  of  solid  angle  centered  at  a .  The  orien¬ 
tational  distribution  of  contacts  A{ot )  is  defined  so  that  A{a.  )duj  is  the  probable 
number  of  contacts  in  this  element  of  solid  angle.  Then  the  average  munber  k  of 
contacts  per  particle,  called  the  coordination  number,  is 


jj  A(a  )du, 


(27) 


where  the  integration  is  over  all  solid  angle.  In  homogeneous  situations,  the  distri¬ 
bution  function  is  independent  of  position  and,  because  a  contact  is  common  to  two 
spheres,  A{—ot )  =  A{oi ). 

When  the  distribution  of  contacts  is  isotropic,  A(a  )  =  fc/47r.  When  the  contact 
distribution  is  not  isotropic,  but  there  axe  orthogonal  planes  of  symmetry  or  an  axis 
of  symmetry,  the  distribution  function  may  be  approximated  by  the  first  two  terms 
in  an  expansion  in  terms  of  completely  symmetric,  traceless  tensors  (  Cowin,  1985; 
Onat  &  Leckie,  1988): 

A(a)  =  +  (28) 


where 


a 


•XJ 


•yl  jj A(,a)(ataj  - 


(29) 


When  the  distribution  exhibits  transverse  isotropy,  a  may  be  expressed  in  terms 
of  the  unit  vector  h  in  the  direction  of  the  axis  of  anisotropy  and  the  strength  £  of 
the  anisotropy  (Jenkins,  1988): 


aij  =  -e{bij  -  Zhihj),  (30) 

where  0  <  £  <  1.  In  this  event,  the  approximation  to  the  distribution  function  is 

(31) 

and  £  is  related  to  the  exact  distribution  function  by 

ei.hh,  -  ift,)  =  jj  Ala  )(a.a,  -  ^6,i)cLj.  (32) 


/l(a)i||;((l-£)  +  3£(/>iOi)"l 


7  Evolution  of  the  Contact  Distribution 

As  grains  slide  at  contacts,  the  orientational  anisotropy  of  the  granular  material 
will  change.  For  example,  an  initially  isotropic  distribution  will,  when  subjected 
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to  triaxial  compression,  develop  transverse  anisotropy  along  the  axis  of  loading. 
The  magnitude  of  the  anisotropy  will  increase  with  increasing  deformation  and  it 
will  be  reinforced  when  the  average  strain  increases  to  the  point  where  contacts  are 
deleted.  It  follows  from  equation  (26)  that  the  sliding  displacement  is  of  the  order  of 
magnitude  of  the  strains  and  from  (25)  that  the  contact  deletion  depends  upon  the 
ratio  of  the  strains.  Consequently,  when  contact  deletion  occurs,  its  contribution 
to  the  fabric  is  much  larger  than  that  of  contact  sliding.  For  this  reason  we  neglect 
the  effect  of  sliding  on  the  fabric  and  focus  on  that  due  to  contact  deletion.  The 
fabric  associated  with  contact  deletion  influences  the  propagation  of  waves  through 
the  aggregate  (Cundall,  Jenkins  and  Ishibashi,  1989). 

There  is  a  relatively  simple  relationship  between  the  coordination  number  ki 
of  a  deformed  aggregate  and  the  coordination  number  k  of  its  initial,  supposedly 
isotropic,  state.  For  example,  in  triaxial  compression,  we  have 

ki=4n  -^sinede,  (33) 

Jo  47r 

where  0i  is  defined  by  (25).  The  integration  gives 

ki  =  k{l  -  cos^i).  (34) 

In  order  to  obtain  an  analogous  equation  of  evolution  for  the  strength  of  the 
anisotropy  e,  we  employ  the  initial  distribution  function  ^4(0: )  =  k/An  in  (32): 

(cos^  ^  —  ^)  sin  6d9  (35) 

o 

In  this  case,the  integration  yields 

£  =  ^  sin^  9i  cos  Oi  =  ^  (1  +  cos  9i )  cos  9i .  (36) 

4  Ki  4 

Once  a  relation  between  the  volume  strain  and  the  shear  strain  is  known,  this 
provides  the  explicit  expression  for  the  evolution  of  the  fabric  with  increasing  shear 
strain. 

We  note  that  in  our  subsequent  calculations  of  the  stress,  we  employ  the  exact 
distribution  function  based  upon  the  deletion  of  contacts  from  the  initial  isotropic 
distribution  function  rather  than  the  approximate  distribution  function  based  on 
the  fabric  tensor  (30). 


HA 

4  ki  Jq 
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8  Stress 

Following  Cauchy,  an  expression  for  the  average  stress  tensor  t  associated  with 
a  homogeneous  deformation  of  the  aggregate  may  be  obtained  by  considering  the 
force  transmitted  over  an  arbitrary  element  of  area  by  pairs  of  contacting  spheres 
whose  line  of  centers  is  cut  by  the  area  element  (Love,  1927,  Note  B).  For  random 
assemblies,  the  stress  is  naturally  expressed  in  terms  of  the  orientational  distribution 
function  and  the  contact  force: 

<0  =  -^  // 

where  the  integration  is  taken  over  the  entire  solid  angle. 

When  the  orientational  distribution  of  contacts  is  initially  isotropic,  the  stress 
may  be  expressed  in  terms  of  the  coordination  number,  the  solid  volume  fraction 
V  =  T:D^n/6,  and  the  components  of  the  contact  force  (6): 

“  4^^  // 

where  the  integrals  are  taken  over  all  contacts. 


9  Pressure 


The  relation  between  the  pressure  p  and  the  strain  for  an  initially  isotropic 
contact  distribution  is,  from  (38), 


where  the  integration  is  over  all  contacts.  We  note  that  the  pressure  does  not 
depend  on  the  tangential  component  of  the  contact  force.  When  (8)  and  (4)  are 
used  in  (39)  and  the  result  integrated  over  ^  from  0  to  27r,  we  obtain 


p  = 


kvM 


fOi 

I  (6  +  ccos^  sinedO. 
Jo 


(40) 


If  there  were  no  shear  strain,  the  volume  strain  Ao  and  the  pressure  p  would 
be  related  by 

fAU 

p.  (41) 


= 


®  kvM 

This  volume  strain  provides  a  natural  scale  for  both  the  volume  strain  and  the 
shear  strain  in  trieixial  compression.  Jenkins,  Cundall,  and  Ishibashi  (1989)  indi¬ 
cate  how  to  obtain  Ao  in  numerical  simulations  and  physical  experiments  from  the 
incremental  form  of  (41)  applied  at  the  initial  isotropic  state  of  the  sample. 
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For  triaxial  compression,  the  integral  in  (40)  gives  (Jenkins,  1988) 


8Ao'^^  =  (56  +  2c)\/6  +  c4-  ^ In  .  (42) 

When  expressed  in  terms  of  A/Ao  and  7/ Aq,.  equation  (42)  provides  an  impHcit 
relation  between  the  normahzed  volume  strain  and  the  normalized  shear  strain  in 
triaxial  compression.  Because  the  pressure  enters  only  through  Ao,  a  single  relation 
between  these  normahzed  variables  applies  no  matter  what  the  variation  of  the 
pressure. 

In  Figure  1  we  plot  A/Ao  versus  7/A0. 


Figure  1:  A/Aq  versus  7/A0. 

With  (42)  providing  the  relationship  between  b  and  c,  we  can  graph  the  vari¬ 
ation  of  6m,  6m  1  and  e  with  shear  strain  in  triaxial  compression.  In  Figure  2 
we  show  the  evolution  of  the  regions  of  slip  and  lost  contact  with  normalized  shear 
strain  7/A0. 

In  Figure  3  we  indicate  how  e  varies  with  7/A0  as  contact  are  deleted.  Values 
of  £  greater  than  one  correspond  to  negative  values  of  the  approximate  distribution 
function  and  indicate  a  breakdown  of  the  approximation. 
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Figiire  2:  6m  and  6i  bounding  the  region  of  slip  from  above  and  9^  bounding  the 
region  of  slip  from  below. 


Figure  3:  e  versus  7/A0. 
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10  Shear  Stress 


The  shear  stress  q, 

9  =  “2(^33  -  <ii),  (43) 

contains  contributions  from  the  normal  and  tangential  components  of  the  contact 
force.  We  first  consider  the  part  q^  contributed  by  the  normal  component. 

Prom  (38), 

9"  =  J f  ^(“3  -  (^) 

Upon  using  (7)  and  (4)  in  (44)  and  carrying  out  the  integral  over  (f>,  we  obtain 

q^  =  f  (b  +  c  cos^  0)^^^  (3  cos^  6  —  1)  sin  6d9.  (45) 

27r  D  Jq 

Upon  carrying  out  the  remaining  integration  and  using  (41),  we  have,  for  triaxial 
compression. 


7^  =  1. 


J  P 
64 


Vb+c 


(362  ^  4^  +  4^2)  _  3  o  2c)  In 


\/c+  y/b+C 

VW\ 


Equation  (46),  used  with  (42)  gives  q^ /p  as  a  function  of  7/A0.  We  note  that 
smooth  particles  can  sustain  a  sheeir  stress  and  that  the  angle  of  internal  friction, 
sin~^{q^ /p)  varies  with  the  shear  strain. 


The  calculation  of  the  part  q^  of  the  shear  stress  associated  with  the  tangential 
component  of  the  contact  force  is  more  complicated  because  its  bilinear  nature 
requires  us  to  employ  our  characterizations  of  the  onset  and  extent  of  contact  sliding. 
When  the  forms  of  the  bilinear  relation  (9)  -  (11)  that  are  appropriate  to  the  elastic 
and  frictional  regions  are  employed  in  (43),  q^  may  be  written,  after  an  integration 


over  (f>,  as 


q^  ^ [  c{b  +  ccos^6)^^^sm^6cos^9(W 

4  A5/2  (2  -  t/)  i  Jo 

—  J  [csin^cos^  — /i(6  + ccos2^)|(6  +  ccos^  ^)^'2  gjji2 
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Upon  carrying  out  the  integrals  using  the  definition  (25)  of  d\,  we  obtain 


P  (1-  i/)  27  f  \/6  +  c,^^2 
48c 


(3b^  +  4bc  +  4c^) 


—  /tsin^xf 


/I  ,  o  M  ( y/c+  y/b  +  c\ 
16c3/a(*’+2‘;)'“(  ^  j 

-  2(6  +  c)(6  +  ccos^  ^jvf )  -  3(6  +  c)^] 

,  ,  .  ^  (6  +  CCOs2^^)l/2  2/>  ^2 

+  H  sm  Om - — - [8(6  +  C  cos'*  $m) 

~  2(6  +  c)(6  +  ccos^  Om)  -  3(6  +  c)^] 


(6  +  CC0S^^A/)^/2  2/>  /l  2/1  V 

- 48c - cos  ^Af  [8  cos^  6Mc{b  +  c  cos^  9  m) 

-  6(6  4-  2c)  (6  +  ccos^  )  +  36(6  +  2c)] 

(6  4  CCOS^  /I  rr,  2  /.  /I  2  1 

4  - COS  ^m(8  COS^  9mC{b  4  C  COS^  6m) 

-  6(b  4  2c) (6  4  ccos^  ^m)  4  36(6  4  2c)] 

62  'I 

-  —  (6  4  2c)  [/ (cos  9m)  -  /(cos  ^m)]  | ,  (48) 


with  9m  given  by  (23)  and  9m  given  by  (24)  when  6  >  0  or  =  arccos  yj-bjc 
when  6  <  0.  Equation  (48)  used  with  (42)  determines  the  part  of  the  shear  stress 
associated  with  the  tangential  component  of  the  contact  force  as  a  function  of  the 
shear  strain.  It  corrects  the  corresponding  exprssion  for  given  by  Jenkins  and 
Strack  (1992). 

In  Figure  4  we  plot  the  and  normalized  by  the  pressure,  versus  the 
normalized  shear  strain. 


Figure  4:  jp  and  versus  7/A0. 

The  relative  magnitudes  of  q^  and  are  consistent  with  Cundall’s  (1988) 
observation  that  in  a  corresponding  numerical  simulation  was  only  about  fifteen 
percent  of  q^ . 

11  Plastic  Strain 

The  average  plastic  strain  associated  with  the  sliding  displacement  has 
the  components 

^  ctj  +  s^Qi)dw  (51) 

where  the  integration  is  over  the  solid  angle  subtended  by  the  region  of  sliding.  We 
note  that  =  0,  so  that  there  is  no  plastic  volume  change.  This  is  not  inconsis¬ 
tent  with  the  relatively  small  inelasticity  in  the  volume  strain  in  a  corresponding 
numerical  simulation  by  Cundall,  Jenkins,  and  Ishibashi  (1989). 

The  plastic  shear  strain  7^  is 

7^  =  -  efi)  =  “  sfai)du;  (52) 

or,  after  an  integration  over  (j), 

7^  =  -  /  [c  sin  9cosd  -  fi{b  +  cos^  9)]  sin^  9  cos  9d9  (53) 
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So, 


7^  =  -  ^c[(2  +  3  sin^  6m)  cos^  Om  —  {2  +  3  sin^  6m)  cos^  6m] 

^\J 

1  3 

-  /i[-(6  +  c)(sin^  6m  -  sin^  6m)  -  — c(sin®  6m  -  sin®  6m)]  ■  (54) 

In  Figure  5  we  plot  7^/Ao  versus  7/A0. 


Figure  5:  7^/Ao  versus  7/A0. 


12  Strain  Hardening 

For  a  given  value  of  the  plastic  strain,  yield  in  the  q-p  plane  occurs  on  the 
family  of  straight  lines 

q/p^hinf^),  (55) 

where  h{‘y^)  is  the  strain  hardening  function.  Because  as  a  consequence  of  the 
kinematic  assumption  (1)  there  is  no  plastic  volume  strain,  curves  of  constant  plastic 
potential  are  straight  lines  parallel  to  the  p  axis.  The  strain  hjirdening  function 
may  be  obtained  numerically  by  first  determining  q/p  from  (46)  and  (48),  then  by 
using  (42)  to  express  q/p  and  7^/Ao  as  functions  of  7/A0  alone,  and,  finally,  by 
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eliminating  the  parameter  7/A0  to  obtain  qfp  in  terms  of  7^/Ao.  In  Figure  6  we 
show  the  normalized  shear  stress  versus  the  normalized  plastic  strain.  This  shows 
the  initial  elastic  response,  yield,  and  the  subsequent  strain  hardening. 


Figure  6:  q/p  versus  7^/Ao- 


Concluding  Remarks 

Although  the  results  obtained  apply  most  naturally  to  situations  in  which  the 
variation  of  the  shear  strain  is  prescribed  and  the  pressure  is  held  fixed,  it  is  im¬ 
portant  to  emphasize  that  they  eire  not  so  restricted.  For  example,  they  also  apply 
when  both  the  sheen:  strain  and  the  pressure  are  varied  arbitrarily.  In  any  case,  the 
relations  obtained  are  universal  in  the  sense  that,  when  the  stresses  are  normalized 
by  p  and  the  strains  are  normalized  by  Aq,  a  single  curve  describes  the  results  of 
all  triaxial  compression  experiments. 

Numerical  simulations  are  in  progress  to  determine  the  ways  in  which  the  kine¬ 
matic  hypothesis  (1),  on  which  the  analysis  is  based,  might  be  improved. 
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Chapter  Nine 


The  Plastic  Anisotropy 
and  Average  Plastic  Spin  for 
Planar  Polycrystalline  Aggregates 

Vincent  C.  Prantil  James  T.  Jenkins  Paul  R.  Dawson 

Abstract 

Significant  changes  in  the  mechanical  state  of  polycrystalline  metals 
accompany  the  large  deformations  associated  with  forming  operations. 

For  example,  crystallographic  texture  is  one  important  source  of  the 
developing  anisotropy  in  the  flow  stress  in  these  materials.  In  order  to 
construct  a  continuum  description  of  such  anisotropy,  we  examine  the 
texturing  of  an  idealized  planar  assembly  of  two-dimensional  grains. 

We  derive  equations  of  evolution  for  the  grain  orientation  and  express 
them  in  terms  of  the  macroscopic  deformation  rate  and  a  single  mi- 
crostructural  parameter.  An  analytic  expression  for  the  plastic  spin 
is  determined  in  terms  of  this  parameter.  We  introduce  a  continuous 
distribution  function  to  describe  the  orientations  of  the  grains  in  the 
aggregate.  We  focus  on  the  second  moment  of  the  distribution  and 
derive  the  differential  equations  governing  its  evolution.  We  derive 
the  equation  for  the  average  plastic  spin  and  express  it  in  terms  of 
the  microstructural  parameter  and  the  second  moment.  This  provides 
a  micromechanical  foundation  for  recent  phenomenological  proposals 
for  the  form  of  the  plastic  spin. 


1  Introduction 

Many  metal  alloys  produced  and  formed  commercially  are  inherently  poly¬ 
crystalline  in  nature.  Each  continuum  point  is  composed  of  many  anisotropic 
grains  or  single  crystals.  When  such  metals  undergo  extensive  inelastic  flow 
and  the  grains  preferentially  reorient,  th^*  response  of  the  polycrystalline  ag¬ 
gregate  can  become  anisotropic.  The  resulting  anisotropy  is  often  referred 
to  as  texture.  The  material  properties  are  defined  by  a  suitable  average  re¬ 
sponse  of  the  collection  of  individual  grains.  The  focus  of  this  work  involves 
using  a  polycrystalline  plasticity  model  to  construct  analytical  expressions 
describing  the  aggregate  texturing  and  the  anisotropic  macroscopic  response 


of  the  material.  In  the  process,  natural  microstructural  variables  are  identi¬ 
fied  that  survive  averaging  and  determine  the  macroscopic  anisotropy.  Here 
we  concentrate  on  obtaining  such  an  expression  for  the  continuum  plastic 
spin. 

We  begin  with  a  brief  description  of  the  geometry  of  crystallographic  slip 
in  a  single  crystal  undergoing  primary-conjugate  multiple  slip.  This  is  fol¬ 
lowed  by  an  outline  of  the  polycrystalline  framework.  A  mean  field  assump¬ 
tion  is  introduced  to  link  the  macroscopic  and  microscopic  length  scales.  This 
is  followed  by  a  relatively  detailed  description  of  the  kinematics  appropriate 
at  the  microscopic  level  of  the  single  crystal.  For  the  case  of  the  kinematically 
determined  multiple  slip,  the  slip  system  shearing  rates  are  determined  and  a 
microstructural  parameter  is  identified  that  enters  naturally  into  expressions 
for  the  plastic  spin. 

The  polycrystalline  aggregate  is  incorporated  by  means  of  an  orientation 
distribution  function  (ODF).  Using  an  expansion  introduced  by  Onat  and 
Leckie  (1988),  a  set  of  even  order  tensors  that  characterize  the  moments  of 
the  distribution  are  identified.  Appealing  to  the  conservation  principle  gov¬ 
erning  the  field  of  grain  orientations  outlined  by  Clement  (1982)  and  Armin- 
jon  (1987),  evolution  equations  for  the  moment  tensors  based  on  analogous 
treatments  by  Advani  and  Tucker  (1987,1990)  are  presented.  In  planar  flows, 
these  can  be  used  to  continually  update  the  ODF.  A  continuum  level  expres¬ 
sion  for  the  plastic  spin  is  derived  by  orientation  averaging.  We  show  how  the 
average  plaistic  spin  depends  on  the  microstructural  parameter  characterizing 
the  single  crystal  anisotropy  and  how  this  spin  is  related  to  the  spin  of  the 
principal  directions  of  the  dev'eloping  anisotropy. 

2  Polycrystalline  framework 

We  investigate  the  material  behavior  of  a  planar  single  crystal  undergoing 
primary- conjugate  double  slip.  The  material  point  is  assumed  to  be  com¬ 
posed  of  a  large  number  of  grains,  each  characterized  by  its  orientation 
relative  to  some  fixed  reference  frame.  The  developing  texture  manifests 
itself  in  a  non-uniform  distribution  of  grain  orientations,  modeled  here  by  a 
continuous  orientation  distribution  function.  We  proceed  from  the  polycrys¬ 
talline  model  and  adopt  a  mean  field  assumption  to  determine  the  flow  at 
the  microscopic  level.  We  use  this  in  conjunction  with  a  multiplicative  de¬ 
composition  of  the  grain  deformation  gradient  to  determine  the  single  crystal 
response.  The  evolving  ODF  is  then  used  to  compute  orientation  averages 
for  the  macroscopic  behavior  of  the  crystalline  aggregate. 

To  construct  the  model,  we  adopt  several  simplif>ing  assumptions  that 
allow  us  to  capture  the  pertinent  features  of  the  anisotropy  at  the  grain 
level  and  understand  how  their  averages  deliver  anisotropic  descriptions  on 
the  macroscopic  level.  We  assume  the  point  is  undergoing  fully  developed 
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plastic  flow  and  neglect  elasticitj"  suppose  the  plastic  flow  is  isochoric;  and 
restrict  our  attention  to  planar  motions. 

2.1  Microstructural  Geometry  of  a  Double  Slip  Sin¬ 
gle  Crystal 

We  consider  a  double  slip  single  crystal  as  described  by  Asaro  (1983).  In  each 
grain  the  deformation  can  be  accommodated  by  slip  on  two  independent  slip 
planes,  each  identified  with  their  unit  normal  vector,  77^“^  Further,  the 
crystallographic  slip  occurs  along  a  unique  slip  direction,  u^°‘\  defined  in 
each  plane.  We,  thus,  characterize  the  crystal  by  the  set  of  time-dependent, 
orthonormal  pairs  imbedded  in  the  lattice.  These  slip  systems 

are  fixed  relative  to  the  lattice  frame  and  to  one  another,  as  illustrated  in 
Figure  1. 

For  convenience,  two  Cartesian  reference  frames  are  employed  to  describe 
individual  grain  orientations  with  respect  to  the  macroscopic  material  point. 
One  is  fixed  in  space  and  characterized  by  the  time-independent  orthonormal 
beisis,  {ci,  62}.  The  second  is  fixed  in  the  lattice  frame  and  identified  with  the 
time-dependent,  orthonormal  basis,  {a,  ax}»  where  a  is  a  vector  in  the  plane 
which  bisects  the  acute  angle  formed  by  the  pair  of  slip  directions.  This  vector 
can  be  used  to  uniquely  describe  the  orientation  of  the  lattice  frame,  and 
therefore  the  crystal,  in  space.  The  initial  texture  of  the  aggregate  is  defined 
by  the  distribution  of  initial  scalar  orientations.  On  for  n  individual  grains. 
To  resolve  the  evolving  texture,  it  is  necessary  to  compute  the  reorientation 
of  individual  grains  with  continuing  deformation. 

2.2  Mean  Field  Assumption 

In  order  to  resolve  the  response  of  the  individual  crystals,  it  is  necessary 
to  understand  the  kinematics  and  constitutive  response  on  the  slip  system 
level.  In  the  context  of  polycrystalline  plasticity  this  is  accomplished  by 
adopting  one  of  several  hypotheses  relating  the  kinematics  on  the  grain  level 
to  those  on  the  macroscopic  level.  Because  we  appeal  to  the  microstruc¬ 
ture  underlying  the  material  point  to  describe  its  anisotropy,  we  must  make 
physically  motivated  assumptions  in  traversing  these  length  scales.  Based  on 
the  classical  work  of  Taylor  (1938),  we  eissume  a  homogeneous  deformation 
throughout  the  aggregate.  The  deformation  gradient  experienced  by  each 
grain  is  equal  to  the  deformation  gradient  experienced  by  the  collection  of 
grains  underlying  the  material  point 

F*=F  =  ||:;i  =  x(X,().  (1) 

Therefore,  each  individual  grain  experiences  the  macroscopic  or  continuum 
deformation  gradient.  This  assumption  ensures  an  entirely  compatible  de- 
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Figure  1:  The  microstructural  unit  model:  a  planar  single  crystal 
oriented  at  angle  6,  the  scalar  measure  of  grain  orientation  corresponding 
to  the  vector  a. 

formation  field.  It  h«is  been  shown  to  be  reasonable  for  moderate  amounts 
of  deformation  in  single  phase  metals  whose  single  crystals  are  characterized 
by  a  number  of  slip  systems  sufficient  to  allow  an  arbitrary  motion. 


3  Kinematics 

At  the  grain  level,  we  adopt  a  multiplicative  decomposition  of  the  deforma¬ 
tion  as  illustrated  in  Figure  2.  The  first  portion,  consists  solely  of  volume 
preserving  plastic  slip,  defining  cin  intermediate  configuration,  B,  in  which 
the  lattice  remains  fixed.  This  crystallographic  slip  takes  place  by  the  rela¬ 
tive  shearing  of  neighboring  close- packed  planes  of  atoms.  In  the  absence  of 
elasticity,  the  remainder  of  the  deformation  consists  of  a  rigid  rotation,  R*, 
of  the  lattice  frame  to  the  current  configuration,  B.  The  total  deformation 
experienced  by  the  grain  is  given  by 

F*  =  R'F**  (2) 


The  rate  form  of  the  multiplicative  decomposition  is  an  additive  decomposi¬ 
tion  of  the  form 

L  =  L*  =  FF-^  =  R’R’'^  +  R*l‘*R*‘^  (3) 


where 


(4) 
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Figure  2:  Multiplicative  decomposition  describing  the  kinematics  of 
crystallographic  slip  in  the  single  crystal. 

is  the  microscopic  velocity  gradient  in  the  intermediate  configuration.  When 
crystallographic  slip  is  the  only  mechanism  for  plastic  flow  in  the  crystal,  the 
velocity  gradient  in  the  current  configuration  can  be  defined  in  terms  of  the 
slip  system  geometry 

lP  =  53  (5) 

0=1.2 

where  7^“^  is  the  rate  of  shearing  along  slip  system  o  and 

0  (6) 

is  the  slip  system  Schmid  orientation  tensor.  Physically,  L**  is  the  velocity 
gradient  accommodated  by  plastic  flow  through  the  fixed  lattice.  In  general, 
the  spin  associated  with  this  pl£istic  flow,  uf^,  will  differ  from  the  overall  spin 
of  material  in  the  grain,  u;.  This  occurs  because  the  crystal  has  available 
to  it  only  a  limited  number  of  deformation  modes  to  accommodate  plastic 
deformation.  The  difference  between  these  rates  of  rotation  then  results  in  a 
subsequent  rigid  rotation  of  the  imbedded  lattice 

w'  =  (7) 

necessary  to  enforce  the  mean  field  assumption  and  maintain  compatibility 
across  grain  boundaries.  The  grain  velocity  gradient  can  also  be  decomposed 
into  its  symmetric  and  skew  portions.  It  follows  directly  from  (1)  that  the 
Taylor  assumption  carries  through  the  rate  form  so  that  from  (3)  the  grain 
rate  of  deformation  and  grain  spin  can  be  written 

d  =  d*'  =  R-(L*’).„,nR*'^  (8) 
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where 


(9) 


'  p 
u;  =  u>  +  u; 

=R’{V'UeuR''^  (10) 

As  a  consequence  of  the  slip  geometry,  the  symmetric  and  skew  portions  of 
the  microscopic  velocity  gradient  can  be  expressed  in  terms  of  the  respective 
portions  of  the  Schmidt  tensors  as 

^  -I- 77*“^  0  1/*“^) 

0=1,2  “ 

=  Yi  (11) 

0  =  1,2 

07^  =  ^  ©  77^*^^  —  7/*“^  ©  I/*"*) 

0=1.2  “ 

=  E  (12) 

0=1,2 

where 

X(“)  =  pb-*)  +  Q(«)  (13) 

Combining  (9)  and  (12),  the  lattice  spin  accompanying  plastic  straining 
can  be  written 

w' =  R-R-^  =  u;  -  E  (14) 

0=1,2 

So,  in  order  to  update  the  grain  orientation,  it  is  necessary  to  resolve  the 
partitioning  of  shearing  rates,  7^“^  among  the  available  slip  systems. 

4  Kinematic  Solution  for  the  Slip  System 

Shear  Distribution 

The  geometry  of  planar  slip  dictates  that  exactly  two  active  slip  systems, 
whose  Schmid  tensors  are  linearly  independent,  are  necessary  to  uniquely 
accommodate  an  arbitrary  deformation.  In  the  plane,  any  additional  slip 
system  will  be  linearly  dependent  upon  this  set.  For  such  redundant  systems, 
there  is  no  unique  kinematic  prescription  of  the  slip  system  shearing  rates. 

For  the  planar  crystal  undergoing  double  slip,  the  available  number  of 
degrees  of  freedom  is  identical  to  the  number  of  independent  components  in 
an  incompressible  but  otherwise  arbitrary  rate  of  deformation.  In  this  spe¬ 
cial  case  there  is  a  one-to-one  correspondence  between  these  two  independent 
components  and  the  shearing  rates  on  the  two  slip  systems.  For  general  iso- 
choric  planar  flows  the  velocity  gradient  has  three  independent  components 
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that  one  can  associate  with  a  rate  of  stretching,  f ,  a  rate  of  shearing,  A,  and 
a  spin,  n. 


L  =  f(t)(e,  ©  ei  -  62  ©  62) 

=  A(0(e|  ©62  +  6206,) 

=  Q(t)(6i  ©  62  -  62  ©  6,).  (15) 

Thus,  there  is  a  unique  relationship  between  {F,  A}  and  in 

(11).  This  limited  number  of  degrees  of  freedom  permits  one  to  kinematically 
determine  the  pair  of  shearing  rates  by  projecting  the  rate  of  deformation 
onto  the  set  of  basis  tensors, 

d0  =  d:  P<^)  =  V“)p(“)  :  p(^)  =  Y  (16) 

Cl  Or 

or 

=  (17) 

0 

where  V~^  always  exists  for  q,  /^  =  1,2  but  becomes  rank  deficient  for  redun¬ 
dant  systems.  Even  though  straightforward,  this  step  has  powerful  implica¬ 
tions  for  determining  the  grain  reorientation  and  updating  the  texture.  In 
this  limit,  the  slip  system  shearing  rates  are  determined  solely  as  a  function  of 
the  macroscopic  rate  of  deformation.  And  this,  in  turn,  directly  determines 
the  plastic  spin  of  material  through  the  lattice  in  terms  of  d 

=  (18) 

01^(3 

From  this  straightforward  relation  the  grain  reorientation  is  cast  in  a  form 
that  explicitly  represents  the  anisotropy  inherent  in  the  single  crystal. 

5  Grain  Orientation  Update 

Recall  that  the  mean  field  assumption  governing  the  spin  of  the  material 
on  the  crystal  level  dictates  the  necessary  grain  reorientation  in  accordance 
with  (14).  In  the  plane,  we  can  recast  this  equation  in  a  form  that  isolates 
the  anisotropy  in  terms  of  a  single  parameter.  For  planar  flows,  the  two- 
dimensional  rotation  tensor  is  simply  a  function  of  the  scalar  6 

cosO  —sinO 
sinO  cosO 
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Combining  equation  (14)  and  the  solution  for  the  shearing  rates  given  by 
(17),  we  obtain  an  evolution  equation  for  the  grain  orientation  in  terms  of 
the  current  orientation  and  the  macroscopic  kinematic  quantities 

d  =  d  —  tja  =  — (19) 

The  plastic  spin  is  determined  explicitly  by  the  relative  slip  system  arrange¬ 
ment  in  the  lattice  frame  according  to  (18) 

—  =  —  ^  =  ASda,  (20) 

a=l,2 

in  which  S  is  the  projection  operator  onto  the  direction  ot^  given  by 

2  =  aj.  ®  ajL  =  I  —  a  0  a  (21) 

and  A  is  a  nondimensional  scalar  microstructural  parameter.  This  parameter 
takes  a  relatively  simple  form  when  the  lattice  reference  frame  is  appropri¬ 
ately  chosen  so  that  ot  bisects  the  acute  angle  subtended  by  the  slip  directions 

A  =  ^  (22) 

^2 

where  Q12  =  QiV  =  Qu  A2  =  ^12^  =  AT  components  of  the 

Schmid  orientation  tensors  relative  to  the  basis  Equation  (19)  can 

be  rewritten  ais 


So 


d  =  d  —  u;a  =  Ada  —  A(a  •  da)a 

=  A[(d(a  0  a)  -  (a  0  a)d)]a.  (23) 


=  A[(q  0  a)d  —  d(a  0  a)].  (24) 


Here,  several  features  of  the  limit  behavior  of  the  microstructural  parameter 
for  single  crystals  give  a  geometric  interpretation  consistent  with  single  slip. 
The  ratio  in  (22)  indicates  A  is  a  scalar  mecisure  of  the  amount  of  material  spin 
relative  to  the  amount  of  shearing  within  a  crystal  subjected  to  an  applied 
deformation  field.  Due  to  the  symmetry  of  the  dual  slip  system  arrangement, 
A  will  vary  continuously  in  a  one-to-one  correspondence  with  the  acute  angle, 
2/3,  between  the  slip  directions  where  1  <A<ooforO</3<7r/4as 
illustrated  in  Figure  3.  Further,  the  variation  of  A  with  ^  indicates  that 
A  never  vanishes.  Thus,  from  (20),  no  crystal  undergoing  double  slip  can 
accommodate  a  prescribed  deformation  without  some  plastic  spin.  Therefore, 
there  can  be  no  inherently  isotropic  limit  for  double  slip  planar  single  crystals, 
t.e.  no  slip  system  arrangement  that  will  allow  the  crystal  to  spin  with  the 
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Figure  3:  Variation  of  the  microstructural  parameter  A  with  slip  system 
separation  20 

macroscopic  spin  of  the  material  point  for  all  arbitrarily  applied  deformation 
fields. 

There  are  then  two  explicit  factors  affecting  the  individual  grain  reorien¬ 
tation.  First,  the  mean  field  assumption  which  directly  affects  the  update 
relation  (19)  through  w.  Second,  the  geometric  characteristic  of  the  grain,  A  , 
which  dictates  the  degree  to  which  the  crystal  will  accommodate  deformation 
through  shearing  and  spin.  With  these  effects  quantified  and  separated,  pop¬ 
ulations  of  grains  can  be  considered  at  the  material  point,  providing  natural 
measures  for  a  continuum  level  description  of  anisotropic  material  proper¬ 
ties. 


6  The  Orientation  State 

Considering  a  distribution  of  grain  orientations  to  be  representative  of  the 
aggregate,  the  ODF  is  the  most  general  continuous  description  of  the  crystal 
orientation  state  underlying  the  continuum  point.  The  orientation  distribu¬ 
tion  /1(a)  is  defined  so  that  A{oc)da  is  the  probability  that  a  grain  with 
director  vector  a  is  in  the  infinitesimal  sector  da  centered  at  a.  Because  the 
subspace  of  possible  orientation  states  is  closed,  the  grain  orientation  can  be 
regarded  as  a  conserved  quantity.  Consequently  its  local  change  is  governed 
by  the  conservation  equation 

/i(a) -H /l(a)div(d)  =  0.  (25) 
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One  efficient  method  for  tracking  the  induced  anisotropy  with  a  limited  set  of 
parameters  has  been  to  expand  the  ODF  over  an  appropriately  chosen  set  of 
basis  functions  (Advani  and  Tucker(1987),  Onat  and  Leckie  (1988)).  What 
results  is  an  infinite  series  whose  tensor  coefficients  represent  progressively 
higher  moments  of  the  distribution.  A  convenient  set  of  basis  functions  is 
provided  by  the  deviatoric  subspace  of  dyadic  products  of  a.  Such  basis 
functions  provide  a  special  symmetric,  traceless  set  of  corresponding  orien¬ 
tation  moment  tensors  which  recover  the  distribution  function.  Expanding 
about  the  isotropic  state  and  retaining  only  the  first  moment,  we  obtain 

Aia)  =  ^  +  -A.,f„{a),  (26) 

2ir  IT 

where 

f,j(a)  =  0,0^  -  (27) 

and 

A,j  =  j>  A{oc)f,j{a)doc.  (28) 

Such  basis  functions  are  appealing  for  several  reasons.  First,  they  are  general, 
depend  on  no  particular  orientation  state  or  frame  of  reference,  and  transform 
as  tensors.  Thus,  for  deformation  paths  leading  to  anisotropy  whose  princi¬ 
pal  directions  may  vary,  one  can  follow  the  axes  of  the  developing  anisotropy 
directly.  Second,  the  lower  moment  tensors  provide  an  approximate  means 
by  which  to  characterize  the  fixed  state  anisotropy  with  a  limited  number 
of  independent  components.  Third,  when  properly  chosen,  they  allow  the 
fully  isotropic  description  to  be  lumped  into  a  single  term  of  the  series.  Fi¬ 
nally,  the  evolution  equations  for  the  moment  tensors  follow  directly  from 
the  conservation  equation  and  their  solution  provides  a  means  to  update  the 
material  state. 


6.1  Orientation  Moment  Evolution  Equations 

One  of  the  advantages  of  employing  the  moment  tensors  to  describe  the  orien¬ 
tation  distribution  is  our  ability  to  evolve  their  components  with  continued 
deformation.  As  outlined  by  Advani  (1987)  there  are  several  methods  by 
which  to  obtain  such  equations.  Here  we  start  with  the  evolution  equation 
for  the  orientation  of  a  single  grain  and  combine  this  with  the  continuity 
equation.  Substituting  (26)  for  the  ODF  and  the  definition  (28)  of  the  orien¬ 
tation  tensor,  an  objective  evolution  equation  for  the  second  order  moment 
tensor  results: 


-I- A(G  +  d)  -  2Aa  :  d,  (29) 
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where 


and 


G  =  dA  +  Ad 


(30) 


a  =  Uijki  e,©ej©eA.©e/  (31) 

a,jki  =  /  yl(a)Q,QjQ/to/da.  (32) 

6.2  Closure  Approximation 

As  a  natural  consequence  of  this  procedure,  however,  the  evolution  equation 
for  the  second  order  tensor  depends  on  the  current  value  of  the  fourth  order 
tensor.  Therefore  a  closure  approximation  is  necessary  to  resolve  a  complete 
set  of  equations.  Various  types  of  closure  approximation  have  been  discussed 
previously  (Advani  and  Tucker  (1990)).  For  illustrative  purposes,  we  adopt 
the  simplest  linear  closure  rule  which  corresponds  to  the  deviatoric  portion 
of  the  fourth  order  tensor  vanishing  (Advani  and  Tucker  (1987,  1990)).  The 
objective  rate  then  becomes 

A  =  A-  u)A  +  Au;  =  ^d.  (33) 

This  is  consistent  with  the  truncated  series  adopted  in  (26).  With  the  rate 
written  in  terms  of  a  microstructural  constitutive  parameter  and  the  kine¬ 
matic  field  variables,  we  now  have  the  means  to  update  the  developing  texture 
with  continued  deformation. 


7  Continuum  Plastic  Spin 

In  discrete  polycrystalline  simulations,  macroscopic  quantities  such  as  the 
continuum  stress  are  evaluated  «is  a  weighted  average  of  the  grain  stress  over 
all  orientations.  The  direct  analog  in  a  continuous  model  is  provided  by  the 
orientation  average.  For  the  plastic  spin  of  the  crystal  this  average  is 

W**  =  ^  A{<x)u^ {oc)da  (34) 

Using  (24)  and  (28) 


WP  =  A(>fd  -  d>l)  (35) 

This  result  has  several  important  features.  It  is  based  on  a  simple  and  definite 
model  of  polycrystalline  slip  and  an  equally  clear  and  definite  characteriza¬ 
tion  of  the  anisotropy  of  the  aggregate.  The  orientation  average  provides 
expression  (35)  directly  and  no  heuristic  averaging  arguments  are  needed. 


The  microstructural  significance  of  the  parameter  A  and  the  second  rank 
tensor  A  are  known;  A  is  determined  purely  by  the  geometry  of  slip  planes 
in  the  crystal  lattice  frame  and  the  second  moment  tensor  can  be  continually 
updated  with  the  evolution  equation  discussed  in  section  6.1.  Questions  that 
remain  concern  the  restriction  to  two  dimensions,  the  accuracy  of  the  ap¬ 
proximation  to  the  distribution  function,  and  the  correctness  of  the  closure 
cissumption  made  when  generating  the  evolution  equation. 

Expressions  for  the  plastic  spin  that  result  from  heuristic  micromechan¬ 
ical  arguments  involving  single  slip  and  dislocation  substructure  also  have 
the  structure  of  (35)  (Dafalias  and  Aifantis  (1986),  Bammann  and  Aifantis 
(1987),  Dafalias  (1990)).  However,  because  the  second  rank  tensor  quan¬ 
tifying  the  internal  structure  is  identified  as  a  back  stress,  it  is  not  clear 
how  these  expressions  for  the  plastic  spin  relate  to  a  plastic  spin  based  on 
polycrystalline  slip.  Equation  (35)  has  the  same  form  as  phenomenologi¬ 
cal  expressions  proposed  earlier  by  Dafalias  (1985)  and  Loret  (1983)  using  a 
second  order  ‘“structure  tensor”  and  an  arbitrary  parameter. 

One  of  the  key  features  of  any  anisotropic  constitutive  model  is  its  ability 
to  determine  the  orientation  of  the  principal  directions  of  the  anisotropy.  This 
provides  one  reason  to  focus  on  the  average  plastic  spin;  it  enters  naturally  in 
expressions  for  the  relative  spin  of  the  a.xes  of  anisotropy  (Zbib  and  Aifantis 
(1988)).  In  fact,  the  orientation  average  of  the  lattice  spin  is  the  difference 
between  the  macroscopic  spin  and  the  average  plastic  spin. 

W  =  <  w'  > 

P 

=  w—  <  w  > 

=  w  -  A(>Id  -  dAl).  (36) 

When  simple  shear  is  initiated  in  an  initially  isotropic  distribution,  the  aver¬ 
age  plastic  spin  vanishes  and  the  eigenvectors  of  the  second  moment  tensor 
begin  to  spin  with  the  macroscopic  spin.  This  is  also  true  in  general  for 
initially  textured  aggregates  when  A  and  d  are  coaxial.  With  continued 
shearing,  the  average  plastic  spin  approaches  the  macroscopic  spin  and  the 
eigenvectors  of  the  moment  tensor  approach  stable  positions  in  orientation 
space  (Prantil  (1991)). 

An  averaging  argument  employed  by  Dafalias  (1985)  has  been  used  by  Ar- 
avas  and  Aifantis  (1991)  to  obtain  an  expression  for  this  relative  spin.  They 
assume  that  the  spin  of  the  eigenvectors  of  the  “structure  tensor”  corresponds 
to  a  weighted  average  of  the  spins  of  “material  fibers”  instantaneously  un¬ 
derlying  the  eigenvectors.  When  the  microstructural  parameter  derived  here 
is  chosen  as  the  scaling  parameter  for  the  spin  of  the  material  fibers,  the 
resulting  expression  for  the  relative  spin  takes  the  form  of  (36). 
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8  Conclusions 


Using  a  relatively  simple  polycrystalline  model,  we  link  the  anisotropy  on  the 
microstructural  level  to  that  on  the  macroscopic  level.  In  particular,  three 
features  of  this  new  approach  are  of  interest.  First,  we  identify  a  microstruc¬ 
tural  parameter  for  the  grain.  It  characterizes  the  single  crystal  anisotropy 
in  terms  of  the  geometry  of  slip  systems  active  in  plastic  flow.  Second,  we 
introduce  the  second  moment  of  the  orientation  distribution  function  as  the 
measure  of  material  anisotropy  and  use  polycrystalline  theory  to  derive  its 
evolution  equation.  Third,  we  calculate  the  average  plastic  spin  of  the  poly¬ 
crystalline  aggregate  and  write  it  in  terms  of  the  moment  tensor,  the  macro¬ 
scopic  flow  and  the  microstructural  parameter.  On  the  continuum  level,  the 
averaged  plastic  spin  provides  the  constitutive  information  necessary  to  track 
the  principal  directions  of  the  developing  anisotropy. 
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